{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><font color=\"darkblue\">Principal Component Analysis</font></h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=https://upload.wikimedia.org/wikipedia/commons/f/f5/GaussianScatterPCA.svg width=250 align=right>\n",
    "\n",
    "### What Features?\n",
    "\n",
    "- High-dimensional data\n",
    "- Data transformations\n",
    "- Interesting directions\n",
    "- Linear combinations\n",
    "- Rotated coordinate system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Location & Dispersion\n",
    "\n",
    "- 1-dimensional\n",
    "\n",
    "> For example, mean and variance\n",
    "\n",
    "- $N$-dimensional\n",
    "\n",
    "> E.g., again mean (duh!) and covariance matrix\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Directions of Maximum Variance\n",
    "\n",
    "- Let $X\\in\\mathbb{R}^N$ be a continuous random variable with 0 mean and covariance matrix $C$. What is the direction of maximum variance?\n",
    "\n",
    "> For any vector $a\\in\\mathbb{R}^N$ \n",
    "\n",
    "> $\\displaystyle \\mathbb{Var}[a^T X] = \\mathbb{E}\\left[(a^T X)(X^T a)\\right] = \\mathbb{E}\\left[a^T(XX^T)\\,a\\right]$\n",
    "\n",
    "> so\n",
    "\n",
    "> $\\displaystyle \\mathbb{Var}[a^T X] = a^T\\,\\mathbb{E}\\!\\left[XX^T\\right]\\,a = a^T C\\,a$\n",
    "\n",
    "> We have to maximize this such that $a^2\\!=\\!1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constrained Optimization\n",
    "\n",
    "- **Lagrange multiplier**: extra term with new parameter $\\lambda$\n",
    "\n",
    "> $\\displaystyle  \\hat{a} = \\arg\\max_{a\\in{}\\mathbb{R}^N} \\left[a^T C\\,a - \\lambda\\,(a^2\\!-\\!1)\\right]$\n",
    "\n",
    "- Partial derivatives vanish at optimum\n",
    "\n",
    "> $\\displaystyle \\frac{\\partial}{\\partial\\lambda} \\rightarrow\\ \\  a^2\\!-\\!1 = 0\\ \\ $  (duh!)\n",
    "\n",
    "> $\\displaystyle \\frac{\\partial}{\\partial a_k} \\rightarrow\\ \\  $?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With indices\n",
    "\n",
    "\n",
    "> $\\displaystyle \\max_{a\\in{}\\mathbb{R}^N}  \\left[ \\sum_{i,j} a_i C_{ij} a_j - \\lambda\\,\\left(\\sum_i a_i^2 - 1\\right) \\right]$\n",
    "\n",
    "- Partial derivatives $\\partial \\big/ \\partial a_k$ vanish at optimum\n",
    "\n",
    "> $\\displaystyle \\sum_{i,j} \\frac{\\partial a_i}{\\partial a_k} C_{ij} a_j + \\sum_{i,j} a_i C_{ij} \\frac{\\partial a_j}{\\partial a_k} - 2\\lambda\\,\\left(\\sum_i a_i \\frac{\\partial a_i}{\\partial a_k}\\right) = $ \n",
    "\n",
    "> $=\\displaystyle \\sum_{i,j} \\delta_{ik} C_{ij} a_j + \\sum_{i,j} a_i C_{ij} \\delta_{jk} - 2\\lambda\\,\\left(\\sum_i a_i \\delta_{ik}\\right) = $ \n",
    "\n",
    "> $=\\displaystyle \\sum_{j} C_{kj} a_j + \\sum_{i} a_i C_{ik}  - 2\\lambda\\,a_k $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And back again...\n",
    "\n",
    "- With vectors and matrices\n",
    "\n",
    "> $\\displaystyle  C \\hat{a} + C^T\\hat{a} - 2\\lambda \\hat{a} = 0$\n",
    "\n",
    "> but $C$ is symmetric \n",
    "\n",
    "> $\\displaystyle  C\\,\\hat{a} = \\lambda\\,\\hat{a} $\n",
    "\n",
    "- Eigenproblem !!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Result\n",
    "\n",
    "- The value of maximum variance is\n",
    "\n",
    "> $\\displaystyle  \\hat{a}^TC\\,\\hat{a} = \\hat{a}^T \\lambda\\,\\hat{a} = \\lambda\\, \\hat{a}^T\\hat{a} = \\lambda$\n",
    "\n",
    "> the largest eigenvalue $\\lambda_1$\n",
    "\n",
    "- The direction of maximum variance is the corresponding eigenvector $a_1$\n",
    "\n",
    "> $\\displaystyle  Ca_1 = \\lambda_1 a_1 $\n",
    "\n",
    "- This is the **1st Principal Component** \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2nd Principal Component\n",
    "\n",
    "- Direction of largest variance uncorrelated to 1st PC\n",
    "\n",
    "> $\\displaystyle  \\hat{a} = \\arg\\max_{a\\in{}\\mathbb{R}^N} \\left[a^T C\\,a - \\lambda\\,(a^2\\!-\\!1) - \\lambda'(a^T C\\,a_1) \\right]$\n",
    "\n",
    "- Partial derivatives vanish at optimum\n",
    "\n",
    "> $\\displaystyle 2C\\,\\hat{a} - 2\\lambda\\,\\hat{a}-\\lambda'Ca_1 = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Result\n",
    "\n",
    "- Multiply by $a_1^T\\cdot$\n",
    "\n",
    "> $\\displaystyle 2a_1^TC\\hat{a} - 2a_1^T\\lambda{}\\hat{a}-a_1^T\\lambda'Ca_1 = 0$\n",
    "\n",
    "> $\\displaystyle 0 - 0 - \\lambda'\\lambda_1 = 0 \\ \\ \\rightarrow\\ \\  \\lambda'=0$\n",
    "\n",
    "- Still just an eigenproblem \n",
    "\n",
    "> $\\displaystyle  C\\,\\hat{a} = \\lambda\\,\\hat{a} $\n",
    "\n",
    "- Solution $\\lambda_2$ and $a_2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA \n",
    "\n",
    "- Spectral decomposition or eigenvalue decomposition or eigendecomposition\n",
    "\n",
    "> Let $\\lambda_1\\geq\\lambda_2\\geq\\dots\\geq\\lambda_N\\geq{}0$ be the eigenvalues of $C$ and ${e}_1,\\dots,{e}_N$ the corresponding eigenvectors\n",
    "\n",
    "> $\\displaystyle  C = \\sum_{k=1}^N\\ \\lambda_k\\left({e}_k\\,{e}_k^T\\right) $\n",
    "\n",
    "> Consider $\\displaystyle C\\,e_l = \\sum_k \\lambda_k\\,e_k\\left(e_k^T e_l\\right) = \\lambda_l\\,e_l$ for any $l$\n",
    "\n",
    "- Matrix form\n",
    "\n",
    "> With diagonal $\\Lambda$ matrix of the eigenvalues and an $E$ matrix of $[{e}_1, \\dots, {e}_N]$\n",
    "\n",
    "> $\\displaystyle  C = E\\ \\Lambda\\ E^T$\n",
    "\n",
    "\n",
    "- The eigenvectors of largest eigenvalues capture the most variance\n",
    "\n",
    "> If keeping only $K<N$ eigenvectors, the best approximation is taking the first $K$ PCs\n",
    "\n",
    "> $\\displaystyle  C \\approx \\sum_{k=1}^K\\ \\lambda_k\\left({e}_k\\,{e}_k^T\\right) =  E_K\\Lambda_KE_K^T$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Combination\n",
    "\n",
    "- If $X$ is a linear combination of the eigenvectors\n",
    "\n",
    "> $\\displaystyle  X = \\sum_k \\boldsymbol{e}_k \\beta_{k} =  E\\,\\boldsymbol\\beta$ \n",
    "\n",
    "> with orthonormal $E=[\\boldsymbol{e}_1, ..., \\boldsymbol{e}_N]$ eigenbasis\n",
    "\n",
    "- We get the (random variable) coefficients\n",
    "\n",
    "> $\\displaystyle  \\boldsymbol{\\beta} = E^T\\,X$ \n",
    "\n",
    "> because $E^T{}E = I$, which is a rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coordinate Transformation\n",
    "\n",
    "- New covariance matrix is diagonal and the elements are the eigenvalues of $C$\n",
    "\n",
    "> If $\\beta=E^T X$ and again assuming $\\mathbb{E}[X]=0$ then\n",
    "\n",
    "> $\\displaystyle \\ \\ \\ \\ \\ \\ \\ \\ \\mathbb{E}[\\beta \\beta^T] = \\mathbb{E}[E^T X\\,X^T E] = E^T C\\,E = \\Lambda$\n",
    "\n",
    "> where \n",
    "\n",
    ">$\\displaystyle \\ \\ \\ \\ \\ \\ \\ \\ \\Lambda =  \\left( \\begin{array}{ccc}\n",
    "{\\lambda_1} & 0 & \\cdots & 0\\\\\n",
    "0 & {\\lambda_2} &   & \\vdots\\\\\n",
    "\\vdots &  & \\ddots & 0 \\\\\n",
    "0 & \\cdots & 0 & {\\lambda_N} \\\\\n",
    "\\end{array} \\right)$\n",
    "\n",
    "<!-- -->\n",
    "\n",
    "> Recall $C\\,\\boldsymbol{e}_l = \\lambda_l\\,\\boldsymbol{e}_l$ for all $l$, so the $(k,l)$ element of the new covariance matrix \n",
    "\n",
    ">$\\displaystyle \\ \\ \\ \\ \\ \\ \\ \\ \\boldsymbol{e}_k^T C\\,\\boldsymbol{e}_l = \\lambda_l\\,\\boldsymbol{e}_k^T  \\boldsymbol{e}_l  = \\lambda_l\\delta_{kl} = \\Lambda_{kl}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigendecomposition \n",
    "\n",
    "- If we multiply with $E$ and $E^T$ from left and right \n",
    "\n",
    "> $ C = E\\,\\Lambda\\,E^T$\n",
    "\n",
    "> or\n",
    "\n",
    ">$\\displaystyle C = \\sum_{k=1}^N\\ \\lambda_k\\left(\\boldsymbol{e}_k\\,\\boldsymbol{e}_k^T\\right) $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dimensionality Reduction\n",
    "\n",
    "\n",
    "- A truncated set of eigenvectors $E_K$ defines a transformation that reduces the dimensionality from $N$ to $K$ while preserving the most variance possible\n",
    "\n",
    "> $\\displaystyle  \\beta_K = E_K^T\\, X $\n",
    "\n",
    "> and\n",
    "\n",
    "> $\\displaystyle  X_K = E_K \\beta_K = E_K E_K^T\\, X = P_K\\,X $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Samples\n",
    "\n",
    "- Set of $N$-vectors arranged in matrix $X=\\left[x_1, x_2, \\dots, x_n \\right]$ with average of 0 <br>\n",
    "<font color=\"red\">*This is NOT the random variable we talked about previously but the data matrix!*</font>\n",
    "\n",
    "> Sample covariance matrix is\n",
    "\n",
    ">$\\displaystyle C = \\frac{1}{n\\!-\\!1}\\ X X^T = \\frac{1}{n\\!-\\!1}\\  \\sum_i x_i x_i^T$\n",
    "\n",
    "- Singular Value Decomposition (SVD)\n",
    "\n",
    ">$\\displaystyle X = U W V^T$\n",
    "\n",
    "> where $U^TU=I$, $W$ is diagonal, and $V^TV=I$\n",
    "\n",
    "- Hence\n",
    "\n",
    ">$\\displaystyle C = \\frac{1}{n\\!-\\!1}\\  UWV^T\\ VWU^T = \\frac{1}{n\\!-\\!1}\\ U W^2 U^T$\n",
    "\n",
    "> So if $C=E\\Lambda E^T$ then $E = U$ and $\\displaystyle \\Lambda = \\frac{1}{n\\!-\\!1}\\  W^2$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm as gaussian\n",
    "from math import pi, cos, sin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzsAAAEICAYAAAB4XuoFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmUXHlWH/jvfbFmLBmRuzJTSm2pkkql2kVVV3f13mxtcHvOweOGsc3mKYPNOcMZPNjgATNgj8Ez2IPdHPfUAIMxbrBnmIYGVwPN0gbc3XSpqlRVKpVU2qXc19j3eL/548bLWDIilSlFZGZEfj/n5FHmey8iX2RV3Pjd33J/YowBERERERFRr7H2+gaIiIiIiIg6gckOERERERH1JCY7RERERETUk5jsEBERERFRT2KyQ0REREREPYnJDhERERER9SQmO9Q1ROTLIvJ39vo+iGhvichPiciv7/ZjiWh/E5HTInJRRJIiYovIT7ThOY+JiBERdzvukXYf/8MRERERUS/4UQB/aox5aq9vhPYPjuxQR7AHhIiIiHbZUQDv7PVN0P7CZIc2EZHbIvIPROQtEYmLyH8UEX/l3H8vItdFZE1EviAiEzWPMyLy90XkGoBrNcf+nohcqwwr/4yInBSRr4hIQkT+k4h4K9cOiMjviciyiKxXvj+8J38EImorEfmHIjJbiQNXReTjIuISkR8XkRuV46+JyJHK9b8gIvcqceI1EfngFs/9vkpMiYnImyLykZpzx0Xkv1Se/0sAhjv/aolot4nInwD4KIDPiEhKRD4nIv+0cu4jIjIjIj8iIksiMi8i31vz2L8iIm9U4s09EfmpPXoZ1AFMdqiV/xbAtwA4DuAJAN8jIh8D8M8r58YB3AHwmw2P+2sAngdwtubYNwN4FsD7oEPMLwP4mwCOADgH4Dsr11kA/m9oz8wUgCyAz7T5dRHRLhOR0wB+CMA3GGPC0JhwG8D/CH3/fxJAP4DvA5CpPOxVAE8BGATwOQD/j9Pp0vDckwD+M4B/Wrn2HwD4LREZqVzyOQCvQZOcnwHw3e1/hUS014wxHwPw5wB+yBgTAlBouOQQgAiASQDfD+AXRWSgci4N4G8DiAL4KwB+UET+2q7cOHUckx1q5V8bY+aMMWsAfhfa6PjvAPyKMeZ1Y0wewI8BeEFEjtU87p8bY9aMMdmaY//CGJMwxrwD4BKAPzTG3DTGxAF8EcDTAGCMWTXG/JYxJmOMSQL4ZwA+3OkXSkQdVwbgA3BWRDzGmNvGmBsA/g6A/9kYc9WoN40xqwBgjPn1SkwoGWN+vvL4002e+28CeMUY84oxxjbGfAnABQCfFJEpAN8A4CeMMXljzJ9B4xkRHTxFAD9tjCkaY14BkEIlphhjvmyMebsSQ94C8Btg+6NnMNmhVhZqvs8ACAGYgI7mAACMMSkAq9BeEse9Js+1WPN9tsnPIQAQkYCI/J8ickdEEgD+DEBURFwP80KIaG8ZY64D+GEAPwVgSUR+szIF9giAG80eU5lK+25lKm0M2iPbbAraUQB/vTKFLVa59kXo6PMEgHVjTLrm+jtNnoOIet+qMaZU87PTtoGIPC8if1qZRh8H8APglNeewWSHdmIO2rAAAIhIEMAQgNmaa8xDPP+PQHtZnjfG9AP4kPOrHuI5iWgfMMZ8zhjzIjSGGAA/B+0cOdl4bWV9zo9Cp8wOGGOiAOJoHgvuAfj3xphozVfQGPOzAOYBDFRilWOqrS+MiHrB5wB8AcARY0wEwGfBtkfPYLJDO/EbAL5XRJ4SER+A/xXAXxpjbrfp+cPQkZ6YiAwC+Cdtel4i2kOVvS8+VokbOej73AbwSwB+RkROiXpCRIagsaAEYBmAW0R+Erqmp5lfB/DtIvLNlYIH/spi5MPGmDvQKW3/i4h4ReRFAN/e4ZdLRN0nDGDNGJMTkecAfNde3xC1D5Md2jZjzB8B+AkAvwXtMT0J4NNt/BX/B4A+ACsAvgbg99v43ES0d3wAfhb63l4AMApd8/cvAfwnAH8IIAHgl6Ex4A+g7//3oNPOcmg+RRbGmHsAPgXgx6HJ0T0A/xOqn2/fBS2asgbtQPm1dr84Iup6fw/AT4tIEsBPQuMS9Qgx5mFmHREREREREe1PHNkhIiIiIqKe1JZkR0R+pbJJ06WaY4Mi8qXKZpJfqqll3vjY765cc01EuP8B0QHGWEJE7cJ4QkRA+0Z2fhW6AWWtfwTgj40xpwD8ceXnOjWL0J8H8ByAf9Iq8BDRgfCrYCwhovb4VTCeEB14bUl2Khu1rTUc/hSAf1f5/t8BaLYT7TcD+FJlE8p1AF/C5sBERAcEYwkRtQvjCREBgLuDzz1mjJmvfL8AYKzJNZOor7Azg/oNKjeIyEsAXgKAYDD47JkzZ9p4q0RU67XXXlsxxozs9X1UMJYQdTHGEyJqlweJJ51MdjYYY4yIPFTZN2PMywBeBoDz58+bCxcutOXeiGgzEdmXu8wzlhB1H8YTImqXB4knnazGtigi4wBQ+XepyTWzAI7U/Hy4coyIyMFYQkTtwnhCdMB0Mtn5AgCngsl3A/idJtf8AYBvEpGByuK/b6ocIyJyMJYQUbswnhAdMO0qPf0bAL4K4LSIzIjI90N3y/5GEbkG4BOVnyEi50XklwDAGLMG4GcAvFr5+unKMSI6gBhLiKhdGE+ICADEmIearronOC+WqLNE5DVjzPm9vo9OYywh6jzGEyJqlweJJ52cxkZERERERLRnmOwQEREREVFPYrJDREREREQ9ickOERERERH1JCY7RERERETUk5jsEBERERFRT2KyQ0REREREPcm91zdARERERES7w7YNZmNZLCXzsAQ4FPHjUL8fIrLXt9YRTHaIiIiIiA4AYwzeuBfDerqwcWw1VcBqqoBzk5E9vLPO4TQ2IiIiIqIDYDmVr0t0HAvxHBK54h7cUecx2SEiIiIiOgDW060TmmZJUC9gskNEREREdAB4XK3X5XhcvZkW9OarIiIiIiKiOuORPlhNWv9ul2A07Nv9G9oFTHaIiIiIiA6APq8L5yYj8LirKYDf48JTR6Jw9+jIDquxEREREREdEKNhP4aDPsSyRVgCRPo8PVt2GmCyQ0RERETU1ZaTeaym83BbgkORPoR8WzfxLUswGPTu0t3tLSY7RERERERdyBiDt2fjWErkN47dXsngzHgYhwcCe3hn+0dvTs4jIiIiIupxy8l8XaLjeG8xiULJ3oM72n86muyIyGkRuVjzlRCRH2645iMiEq+55ic7eU9E1H0YS4ioXRhPqJcsJTcnOgBg28Baj+6bs1MdncZmjLkK4CkAEBEXgFkAn29y6Z8bY76tk/dCRN2LsYSI2oXxhHqJJYKybbCYyCGe1Q1DI30ejPb7YPVuzYEd2c1pbB8HcMMYc2cXfycR9R7GEiJqF8YT6mqjYS9uLqewlMwjX7KRL9lYSuZxZzVzYAoQ3M9uJjufBvAbLc69ICJvisgXReSxXbwnIuo+jCVE1C6MJ9TVDAQhf/1ELZclGAn7sMppbAB2qRqbiHgB/FUAP9bk9OsAjhpjUiLySQC/DeBUk+d4CcBLADA1NdXBuyWi/YqxhIjahfGEekEyV8R4pA+DQS8S2RIsSxDpc8NtWUjmihjr9+/1Le653RrZ+VYArxtjFhtPGGMSxphU5ftXAHhEZLjJdS8bY84bY86PjIx0/o6JaD9iLCGidmE8oa7n97gAAD63CyNhH4aCXrgta+MY7V6y851oMUwsIoeksm2riDxXuafVXbovIuoujCVE1C6MJ9T1xvr98Hk2N+c9bguHIluP6hhjkM6XkC+VO3V7+0LHp7GJSBDANwL4uzXHfgAAjDGfBfAdAH5QREoAsgA+bYwxnb4vIuoujCVE1C6MJ9QrXJbgmakBvDufQCxTqcYW8ODMoTA8rtZjGvPxLK4vpZAv2hABRsI+PDrev+VjulXHkx1jTBrAUMOxz9Z8/xkAn+n0fRBRd2MsIaJ2YTyhXhL0uXH+2CByRR2hcaa2tbKeLuCd2cTGz8YAS4k8ynYcT08NdPRe90LvpW9ERERERAeM3+O6b6IDADPr2abHV1MFZAqldt/WnmOyQ0RERER0QOS2WKOTK9q7eCe7g8kOEREREdEBEenzND1uWUDItyu70uwqJjtERERERAfE1GAAHvfmFGBqMAhvk+PdrvfSNyIiIiIiasrvceEbjg3g1koasUwRXreFiWgfJqN9e31rHcFkh4iIiIhoH8kWypiPZ1GyDaIBD0ZCPlS2fmqLgNeNxyYibXu+/YzJDhERERHRPrGUzOHSbBx2pVbA3VVgKOTFk4ejsKz2JTwHRe9NzCMiIiIi6kK2bfDufHIj0XGspgqYizcvGU1bY7JDRERERLQPrGcKKJaal39eSuZ3+W56A5MdIiIiIqJ9wNpiXc5W56g1JjtERERERPtANOCBz9O8eT7W79vlu+kNTHaIiIiIiPYBEcG5iQjcrvpRnPGoH+OR3iwN3WmsxkZEREREtE8MBL14cXoYi8k8iiUbgyEv+v2evb6trsVkh4iIiIhoH3G7rJ7d5HO3cRobERERERH1JCY7REREREQPyBiDeKaIZK6417dCTXAaGxERERHRA1hO5nFlIYF8UffGCfrcODfZj/B91tispvKIZYvwuS2M9fvhcXH8oVOY7BARERER7VC2UMbbszHYNXuApvMlXLwXwwdODsOyNu+LY9sGb9yLYT1d2Dh2fSmFp6cGEOljEYJOYBpJRERERLRDs7FsXaLjyBdtLKfyTR9zdy1Tl+gAQKls8NZMDKupPKfCdUDHR3ZE5DaAJIAygJIx5nzDeQHwCwA+CSAD4HuMMa93+r6IqLswlhBRuzCeUDsUSk0ynfucW0puToIWElksJfOIZ4rwe1yIBjw4NxmB3+Nq270eZLs1je2jxpiVFue+FcCpytfzAP5t5V8iokaMJUTULown9FAGgh7MxbJNz0UDzaekGWPqfo5lC1iI5yvnKscyRbwzl8CzRwdQLNsolQ36vEx8HtR+WLPzKQC/ZvS//tdEJCoi48aY+b2+MSLqKowlRNQujCd0X2NhP+71ZZHI1k89OxTxtyxQMNrvRzKX2vh5tTKlze+x6hKa5WQOr95eRSJbgjFAwOvC9FgIo2F/B15Jb9uNNTsGwB+KyGsi8lKT85MA7tX8PFM5VkdEXhKRCyJyYXl5uUO3SkT7GGMJEbUL4wntSDxbxL21DJYSOdi2DsFYluCZqShOjoYQCXgwEPTg8GAfjDH42s1VvD0TR7whEZoaDGAgWE2ESmUDtyU4MhCou+7OagZz67mN0Z5Modz0+ej+dmNk50VjzKyIjAL4kohcMcb82U6fxBjzMoCXAeD8+fPmPpcTUe9hLCGidmE8oW2xbYNLc3EsJaprbfweF56eiiLoc8PtsnB8OIjjw0GspQu4eG99o2hBKlfCciqHJw9HMRTyAQBcluCZqQGspAqIZ4vwuAXZQhluqzr+kCuWkS2WEWiYumYMcG8tg8hkpPMvvId0fGTHGDNb+XcJwOcBPNdwySyAIzU/H64cIyLawFhCRO3CeELbNbOerUt0AE1G3plLbLr2xnJqU3U22wZuLKfrjokIRsI+TI+G8PzxIYR89VPeCiUb45E+aJ2MevlS+QFfycHV0WRHRIIiEna+B/BNAC41XPYFAH9b1PsAxDknlohqMZYQUbswntBOLCRyTY8nskVkCqWNn40xiGeaTzFLZIsbU98auS1NfNKFElZSebhdghemhzAc8ja9/n6bldJmnZ7GNgbg85XM1A3gc8aY3xeRHwAAY8xnAbwCLe14HVre8Xs7fE9E1H0YS4ioXRhPaFtyxTLmYlmk8iX0+90IeOubzbX5i4jA67aalpz2uK2mG4yWbYMLd9aRypUQ9LoR9LpRKhukciUcHgjg3lpm0/NMDQY2PQ9traPJjjHmJoAnmxz/bM33BsDf7+R9EFF3YywhonZhPKH7SedLeP3uOt6ZiyOeKaJkG/jcLoyEfJgc6AOg1dFCvvpm9OGBPtxsmLIGAJPRvqa/Zy6WRSpX2nT87loG7z85hIDXhdlYFsWyjYGAFydGgtx75wHsh9LTRERERER7LpEr4i9vruHt2RhsGxDR8tADfV4sp/Lo73MjGvDizHj/psceHw6iULYxF8vCtgHLAsYjfTg5Emz6u2Itpr0ZA8SyRRwZDOAIR3IeGpMdIiIiIup5uWIZ1xZTWE7lIBCM9vtwajQMr7u6hP3GUgrr6cJGoQFLLEz0+1GyDQaCHgR8LrxwcmjTCEsyV8RcLIdS2eD4UAiRPg9CfnfdczfyuDdPbXN4XbuxO8zBwGSHiIiIiHpaqWzjwu115IpONTOD+VgOyVwJzx8f3Kh8tpLKYymZw3w8C0uAkN8DtyVYSxcgInBbFm4up3FqLARPJSGZj2dxeS6xsScOAHjdFqbHguj3ezdNd3NMRPswu56texwA9HldGAw2L1BAO8dkh4iIiIh62kIiV5PoVKVyJaykChgJ+1Aq27i1kkEmX0K+ZMMYIJHLIl8qI+TzwBKg3+/GXCyLbLGMZ48OoGwbXF1IbiQsBgb31jJYz+gmpGP9fgyGvHh8MrKRHDn6/R48Ot6P9xaTKJX1CYI+Nx4/HGladpoeDJMdIiIiIupp6Xzr/WnS+RJGwj7MxrIIel1IZC0MB31YSeeRLegGnyMhP0b6fRuln9fTuimobZuNRAUAlpJ5rKV1LU4iV8RYvx9rqQKuLiRxrslmoBPRPoz1+xHPFuGyBJE+lpZuNyY7RERERNTTAt7WVcwCPj0XyxQxGvahULYhAvg9Fu6UbQR9bpyb7MdktL5YQLZQRl/D866lChvfu2pGZ5aSOZTtfrialKB2WcJpax3EZIeIiIiIetp4xI/bq2nki/X74AR9boyEfAAAj8uCiODIQACH+v3IFcsYDfuxlikg0rc5GQn6XAj7PQj63EjntYR0qWbznYGaBMa2gZJtw2WxdPRuY7JDRERERD0hXyrj5nIay8k8RICxfj+ODwfhcVl49ugAri2msJLSc6NhP06NhTbWx0wO9GEulgWgiY/HZcHvcSFfLm8qMjAcrk5pe/xwBBfvxpArlhH2uxHPFjEc8mEgUE12Qn43fG4mOnuByQ4RERERdT3bNnjtzjoyNetz7q5mEM8Wcf7oAAJeN548EoVtG4hgUxGASJ8HZyfqCwYMBr04f2wAc7EcVtN5uCwLExE/ToyENh4X8rnxgekhrKYLODYcwI3lFFxSLUZgWcCp0RBobzDZISIiIqKut5jM1SU6jnimiNV0AcOV6WpWk3UzDqdgQCJbhNslG6M3I2H/lr9bRDAc8mE45MPhgQBm1jNI5Ero87hwZDDQsvw0dR7/8kRERETU9ZK50pbnnGTnflyWNKy3Mbi3nsFCPAfbACNhL44OBTeVknb4PS5Mj4Z3dvPUMUx2iIiIiKjr9Xk2r4nJl8pYTuVRNgapXBFTg0FEAjsr7/z2bBzLyfzGz+m87s3z3LHBLUeJaH9gskNEREREXe9QxI9bK2kUSlpxLVcq49piCslcEbF0AW/PxBHwuvDR0yN4dGLznjfNxDPFukTHkcqVsJDIYSLa19bXQO3HZIeIiIiIup7HZeGZowO4upDAerqIxXgO2UIJHpeFSv6DTKGMP76yhLFI37b2tolni1ue2yrZWUsXcHM5hXi2CJ/bhcMDfTg6FNhUGIE6q/lkQyIiIiKiLhPyufHs0UF86JERHB0KoM/r3rS2Jle08d5iclvP5/e0bir7m0ybc8QzRVy8t45YpghjgFyxjOtLKVxfSm3vhVDbMNkhIiIiop7idVtwWRbKNZt8OiwLyBc3V21rZjjka5rUuCzBeKR1hbY7a2nY9ubjM+tZFMtNTlDHMNkhIiIiop5zfDgIV5MCApE+D8J92ytSYFmCp6ei6K+5PuB14akj0S1HdlL55pXhyrZBbpuJFrUH1+wQERERUc85NhzE44cjeGsmtjHKEunz4MhAAMeGgtt+nqDPjeeODyJbKMM2BsFt7JkT8Lqb7vljWVtPf6P269jIjogcEZE/FZHLIvKOiPwPTa75iIjEReRi5esnO3U/RNS9GE+IqF0YT7pfOl9CusXISaOPnxnFtz8xgbMT/Tg73o/HD0fw1FR0W8UJGvV5XdtKdADg6GAAzeoQTEYDLffnoc7o5MhOCcCPGGNeF5EwgNdE5EvGmMsN1/25MebbOngfRNT9GE+IqF0YT7pUIlfE5bkEUpXNQ4M+N86O92+5b46I4NRYGKfGwrBts+19cdbSBdxdyyBTKCHs8+DocAD9/u3vzzMQ9OKJw1HcWE4hlSvB7RIcHgjg5Mj2R5SoPTqW7Bhj5gHMV75Pisi7ACYBNAYTIqItMZ4QUbswnnSnUtnGG3djKJaqi/vT+RLeuLeOD0wPb2u0ZLuJzlIih7dn4zCV2gaZfBlLySzOTUYwEvJv+3lGwj6MhH0o2waWgCWn98iujKOJyDEATwP4yyanXxCRN0XkiyLy2BbP8ZKIXBCRC8vLyx26UyLa7x42njCWEJGD8aR7LCRydYmOo1Q2WIjn2vq7ri+lNhIdAFhO5nFpNoHfvjiHP7u2XDm/ucpbKy5LmOjsoY4nOyISAvBbAH7YGJNoOP06gKPGmCcB/BsAv93qeYwxLxtjzhtjzo+MjHTuholo32pHPGEsISKA8aTb5JskOts5t1OFko1MoVpYYC1dwGwsi2LZIJMvo1Q2uL2Sxs2VdNt+J3VWR5MdEfFAA8l/MMb8f43njTEJY0yq8v0rADwiMtzJeyKi7sR4QkTtwnjSfSJblIpudc4Yg/V0AavpPOwm++04bFuvi2UKcAngclVHYZaT+Y3v3TXT1+6tZbZ8Tto/OrZmR3S87pcBvGuM+ZctrjkEYNEYY0TkOWjytdqpeyKi7sR4QkTtwnjSnYZDPgwEvVhPF+qORwMeDIc2V1abWcvgT64ubSQrwyEvPnp6FFMNJadXU3lcnk8gX9TRIZ/HQsjnRjxTBAAUajYAHQn7Nr4vlQ2Ktg2fxTLS+10nq7F9AMDfAvC2iFysHPtxAFMAYIz5LIDvAPCDIlICkAXwabOTSZBEdFAwnhBRuzCedKmnjkRxdy2DxYSu0RkN+3B0KLhpPcx6uoAvXlpAPFvcOLacLOCL7yzgO545jKGQJi35UhlvzcRRrhmhyRdtFEs2RsI+rKbz8Hss5Io2hkNeDIeqyY7f44KXJaS7Qiersf0FgC1XYxljPgPgM526ByLqDYwnRNQujCfdy2UJjg8HcXx46/LNN5ZTdYmOI5kt4fpyaiPZWYzn6xIdh22AwaAXj47348RICO8tJGE1JFTHRzYnWbQ/dXJkh3pAJgMsLuqOv+PjgHebe3BlMoAI0NfX2fsjOmiKRWBhQf8dHQVCob2+o4dXLAL5PBAIaKwhInoYidzmRAcAbGNjZi2DR0bDiAY8dVPUGhXKNrxuC8eGggj73Li9mkYyV0LA68bUYACHIv5O3T61GZMdaun6deDKFWyUX7x0CXjqKWBysvVjYjHgzTeBRAIolfTr5ElgagoY5tJOooeysgK8+qq+rwDgnXeA48eBc+f29r4eVKkEvP02MDurccbnA86c0XhRyxhN8BIJIBgEJiaYFBFRa4f6/bBER2gciVwRsUwBQ0EfXruzjoDXhcmB1j2yg4Fq7+5QyLcxGkTdh8kONRWPA+++W3/MtoGLF4GRkc0jPPk8kEwCX/86UC7r9++9p9+/+y7wxBPaQDl/no0Uogdh28Brr1UTHcetW/qeHBvbm/ty2DaQSmnC4ttmm+CttzTRceTz2lni81VfT6EAfPWrmug4rl4FXnhBR4KIiBpNRPsQ6fNgOZmH22UhVyxjNVXAUNCzMSKTyBWRX7YxFPJiNVVf9GC0X4shUG9gskNN1TZAatk2MD8PHD2qPxeL2jhZWNDjc3Oa1CwuaqLjXLO+DrhcwJ072hNNRDuzvKwN/2ZmZ/c22bl9WxOQQkGnr46PA08+Cbi3+ITJ5zVeNHPrVvX1XLlSn+gAOk320iXguefacvtE1CNs2+DdhQQW4jkciviRzBWRzJdhCTA9GsIjoyFYluDWShqJXBHGAI9N9OPIYGBjSttYvx8TnKLWU5jsUFP2Fvtz1Z57801NcgBt6JRK1RGdYM36QaeRNj/PZIfoQWxVB2qr92unLSzoVDSHMZrEGKMjua3kcq1fUzZb/b5VQrS0pHHGxaqvRFRxezWN+ZhWavO5XTg3GQUApPMlBH3a5H1vMVm3aWg6X8JyMo/nTwwi4GWzuBdxQhE1dehQ8+Mi1XO5nDZ0HE5yY4yu3alVe46Idm54uPVISav36264fbv58YUFjRGthEKtX8/AwEPfFhEdQLOxak9JybYxF8viykICl+cTWIhnkcwV6xIdS4Cgz42ybTC7nm32lNQDmOxQU8PDmxcJA7p42Kmw1tgzOzCgc+iLRWB1VaeyZbNAfz8Qiehxy6r2yJbLm5+fiJpzu4HHH9cOh1pjY1sXDem0VgmNMTpVrRWXCzh1Sqe43rypU9ficX2d09PV68bHmz9+ZGR7ozrOdLnl5ft3tpTL7JAh6malsr6BbdvgxlIaS8k8ckUbbkuwnini1kq67vrRfj88lb1yskU2SnoVx+uopSef1EbUwoImKZOTmrQ4wmFtmDgLpkW08RGP679OUhMOA7//+1pJ6tQpIBrVaW3Hj+v6njNn9n5xNVE3OHxY3z8zM9XS06OjmxOg3TQwoAVJGnk89y+LnUzq1507GisyGeD06frHnTkDrK1p8QNHX199BbqlJZ0+u76usWViQuPX7Kwed6b59fXpOp/+/vr7WF7WtUGxmMa0qSng0UdZTIWo2wwEvVhJ5rGeKdQlL0GfG9OjIczFsohniwj63BgMeNHf59m4Juz3NHtK6gFMdmiTQkETlr4+HeFpVTLa5QIeeQS4fFl/XlvThkWhoEmNy6W9tV/+sn4fCmnZ3EJBGyKlko4Evfoq8P73A4ODu/YSifaN9XVt6A8Obq9xHQppArBfTE/rWrxiw7YWTgxoZWUF+MpXtDPFmeaaTgOvv66vLxzWYz4f8OGh2ApfAAAgAElEQVQP6+9wSk9PTupzF4ta7fGNNzRJuXFDj128CHzta1o1snaEOpvVePOxj1UTxHhcq0g6CVGppCNNhQLw9NPt+RsR0e44ORLEeqZQP1XNAsajfbBEcHgggMNNpsn6PBYmo9wYsFcx2aE6V65og8H54B8eBp59tvVmoidPalJ0+7Y+VgQ4ckSvL5W0B3p2VhtymYz2zoZC+jvCYW28RKP6M5MdOkgSCeDCBW3gAzoS8vjjezsl7UEEg8CLL+q+XGtrgN8PHDumoytbee+9+jV/jpkZjSePP1495ows1/5t7t3TwgiXLunf8u5dYGioOip09arGr3C4fg1QJqPTbJ1OnJs3mxd4mJ3V0R0/izIRdY2w34Pnjw/Ctg3y5TJ8LheGw170earN3ccnI8gUy5iPZVGyDYZCXpwYDsHr3v5Qrm0bzMayWEzkYACMhHw4MhiAy9rDYXZqickObZiZAa5dqz+2sqK9pFuVeJ2Y0K/Z2foSsUtL+rMxmvh4PPpvOq0NiEKhum6ndopKI6fH2MMRZuoRtg385V/Wr3cpFnWEIhzePM1qvwuFdMPhnWgsYlJrbW3rx6bTWgnSGE1eslldm+NUe3S7q/FlZWVzwYPavYpaxR5jqrFqv7Ft/f/F693bKYxEu6FUtnFrJY3FhC4CHO334fhwcGOtTaOA140PTA/DEkHZrl+E1+d1YSTsg2UJjg8Hmz5+O96ejWM5WV2UGM8UsZrO45mpAQjflPsOkx3acOdO8+OLi9oou9+H/qFDukmgI5vVD2WPp/5DuVjU742pNuqaNe5SKe25XVnRn0dHtbeXGwlSt1taar6w3xgdsXjssd2/p902NlaNA7VE7r+Gb3a2+ri+vmonizG6BmhgQGOKy7W5EIrLpSNAjv7+5omXSH35/P3AGB1Bv327Og349Gldy0XUi4wxeONeDPFMdZ7s3dUM1tIFPHdsEFaLkRS/x4WnjkTx7kICmbwGgWjAg7MT/S0fs12xTKEu0XGsp4tYTuUxGt6HPSQHHJdf0oatKie12syw1qOP6oeuk9SIaEGDYFATnFhM/y0W9fiJE5oIWZZOh6tVKumu6U6iA2gD8atf3ds9RYjaYav303bea73g1Cmd7lbbCepUaDt2bOvH1o7MHDqkCU+xqIUGbt3SabH9/TriFI3WP/bMmfpR4hMnmq8tOnJk/43qXL6s0wWd15/J6Gjg4uLe3hdRp6ykCnWJjiOVK2GpScJRayDoxftPDuOFk0P4wPQwzh9rzz46603uxxHb4hztHY7s0IbhYR1NaRyB9fk2V1WKx7UHNRisThEZGtJCA++9p+d9Pj139ar2Yns8mqj4fDot7vBhPX/mzOYGyeysPiaR0MZfMKgNmrU1nVY3Pd19U33oYDGm9RSj2pGFnZzrtFRKOyX8/taFSXaiXNbOjGZ/h/5+4EMf0vd+PK7XDAxoUYD7JRmjo5rQALrWb3xcCxUUi/r3i0Y11oyNaXGD1VUdXT5yZPPawHAYeOEFHTFZXdX4dPSoJl37SToNvPaafh+J1BezuHFjexUtczktiOH3cy8j6g7JXOvkIZkr4lDk/j0Szmai7bLV2h7fDtb90O5hskMbFY1u3NBpY8GgNgqc/XTOnq1+sJbL+oFb25MYjWoD4c03tbHU368JzAc+APzO7+gHczyuvZCWpT2xpZI2PoaGqh+6tq3JjGVpo+PSJX0MoA3HbFbvaXVVk6GREd2hvdXGhER7IZ3WHvjFxerC+rNn60cTgkFdW3LrVv1jo9G9mZJkjHYizMxUj4XDwPPPV+PATiwv65qkmzf1b/DYY5rYNI6gTE9rorKwoMnOxMT2RlOGh/XvOjurPxeL+jtCIZ3aVS5rHBoe1r/niRPNnyceryY5Xq9OCZuebv86mJ1W3Gt086aO4DhrKp29iJwOn3S69WMdly/r8zjT/yIR7XTab6NXRLX6vK1LOvo929hoqwPGwj5cc8nGnj4OlyXbSr5o97GZSPj61zXJcLu1wbCwoI2I979fezdre0KvXt08ZeK//lf9EHV6FtfWdBO/7/s+TXoWFvTDWEQbQSI6ZS6V0g/vVEobLm+9VZ3Cc/Gijuo4H8Rra9ogGRnRcteAPtelSztfGE3UKcWillN21uOUy1olLJEAPvjB+mvPndNEf3ZWk/+xMZ2+tRd7u9y6VZ/oADpye/GijnrsRDwOfP7zuvbIsbioMeG7vmtzIhEMbp7Guh3PPKPJ0fy8Pv/w8OYRYkBjT7NqkqmUxi5nTU82q4lPLldfCe5hxOPaOfQwFffW14F33tH47KxBKpV0OtuTT+qx2v3Pmrl3rzoSVntvr7+ucZ5ovxoN+3Hdk0K+WD9/3eO29iyxcLssPD01gHdm4xslrv0eF85O9MPn3psEjLbGZOeAW12tr3xUuy/F0FB9opPPa+9isaiNCmefi7ff1sacz6e9wS6XJi1f/rKOEDmNm5mZakPOWauzuqqJUjBY3RyxUNAGx8qKNmYsSz+YAT1X26CZnQWeeIKb/9H+cO9e88IDsZgm5yMj9ccbyyl3Qqmk691sW99jzRr+tYlJrZWV+xcnyeerJaTHxrRh3pg4AdpRcufO/dfj7MShQ/pl25pMNbKs1kUGbt7cXLwA0OT0kUc0ntUqFvV35POaWN2vVL5ta0dSs4p7/f3VfYTux/lbWpaOgjk/l0oaF4eHdZRnK3fvNj++uqqJ2H4rxEDdxbYNYtkiXCKIBNpbNtVlCZ6ZGsDluQTurmVgYDAZDeDcZH/Lamy7IdLnwfunh5HIFbXYkt/NKmz7GJOdA26rks+1u6LfuqXTIJzdyF0u/YCdm6t+kC4uaoNufFynk8zO6uZ9X/+6Tp1IJLSX0hidwnPpkjYcVla0YTE2ptNIbLs63c3n0+stS5Ocycn6nmHbrq4LINprte+ZZucak51OW1rS3nunfLszpawx4WjW6HfUFgNodO+ejsg6RUNENKFprLAG6LG7d9uX7Ni2PqfLpSNDCwubi5dMTWlyl83qCFC5rHHF2eOr1fMmk/XJzuqqxjHnb3H1qsa5Z59tPeXNqWLZyKm4d/bs9l5n7d9/YkI7ihYWtOMnEgHe9777J15b/Tfc6hzR/SzEc7iykNiY0hXwunDucAT9/vYlPcWyjWyxjEBlSlu2WMJ6uohwG3/Hg2rn66TO6XiyIyLfAuAXALgA/JIx5mcbzvsA/BqAZwGsAvgbxpjbnb4vUo2FB2o5PY/xuCYmgM6Hv31bG0+zs9pDKaJTLJwysgsLOkc+ENAP/FhMp4e8+65+UI+NaeLkduv35bIez2S01/L48WpDY3y82jhZX9/8oR6Ncv+dg6IbYslWvfXb7clvl2JRp1DVNmZtW0diBwf1/MxMNUFoVlAhGGwdI7LZ6l43DmO086JQ2DyCJNKev0E+r/HISW5GRnRK4PPPa5xZX9fffeyYjtDcu1d/n1eu6PTcYFCvbdRYctoYHY1pTArm5zV5O3q0+X22q+Le6Gj9SNnIiH6JAJ/4xPbXODVL7pzR+IOoG+LJfpfKl/DOXLwuBmQKZbx5L4YPnBx+6BLPgI4avTUTR6Fkb4zk2Dbw3mISkT7PA40k2bbBQiKHtXQBlggmon5EAy12Tqee0NH+cBFxAfhFAN8K4CyA7xSRxv6s7wewboyZBvCvAPxcJ++J6tUWCKhVO53N+aBNJjVxSST0+7t3dR740JB+YGazOiUil6tuxnfvnvY++nw6auP1ahKUyWgSlUrph3nt2hxAGxAi1QW4k5P6uPHx6j26XNvvHaXu1i2x5MiRzdOfAE3Kd3tUZ2Ghda/9n/+5ri26e1ff37GY/lvbaLEsTSJamZurXm/b1aTp1CmNBU65+URCzw0Pty4UsF3GAF/7mv5u5/ctL+triUSAF18Evu3bgG/+Zo03hYKOPDWONF27pv9Nmo0Ij4/XF2VYX9fX0+pv0MpW1ex2UuluYqL5/zunT2+/uMD09Ob9yUR0lO8gjop3SzzZ7+Zj2aajuPmijZX01mWht5IrlrGUyCGeKWIlnUeh1Hy/iflEizfmFmxb9+25PJfAQjyHuVgWF26v487qNqp8UNfq9MjOcwCuG2NuAoCI/CaATwG4XHPNpwD8VOX7/xfAZ0REjGn2FqJOeP55HWmZna1Oj5iY0Ln3jzxSneJy546OxkxNaaLi9N6Gw9pbGovpKEskooUJnBK68bg+h9eriU0qpT+7XPp8Tzyh0+Rqp9RFo8C3fqsmO7lc9efZWW08BYPac7vVyBT1lK6IJR6PViF85x2dQiZSrcZWa2lJOwqc6oXT0+0vOd0q0XHW2NSOSHg82sgfG9P3pN+v57d6f9m2vjfv3q2uqYtG9XFPPql7YqVS+j63beBTn7r/Qvr7WVpqPkKRz2vHyokT9aNTjVPbSiWNWz6fdrh8wzfoiHMiofd55Mj9O1CM0ee5XxVIJ0bdvl1/fHBQ4+t2iWjVtLk5ff0ul1aY28n/Lz6fFsi4c0c7lJz/vs0KOhwQXRFP9rtCufWmd60SlPu5spDA7Ho1iSqUyxBI0/U5jdXQtmMhkcN6evPQ6o3lFMYjfVuWlabu1elkZxJA7dLXGQDPt7rGGFMSkTiAIQArtReJyEsAXgKAKWfIgdrC49HGidutH+K1BQVmZ/UDO5erloF2u/WD1u3WxsSdO9qgGxrSER/L0p7Qr3xFe4aNqU5vicWqU2PKZX3Ot9/W5xXRD+DRUX2+ZiV4ubfOgdU1sSQY1Map82HdODVsbq66Xwqg763lZe102Gr0xxgdMXW5tlcOenS0+fFYrPn7yO3W533iifs/N6Dvd2dvG8f6uq5vOXtWR1jicU02IhE97sSXGzc0thijSdb0dP10VKeDpHH0ZasSy1utP3TWCy0tVdcAWpbe5+iovgaXq/kox8CAJge5nMa7+Xm93uPRPXy28vjjmtzMzurrGRvTGLfT0RTL0nj4MGXJvd79t3fQHuqaeLKfDQS8mI81WZgGYDC482lhM+sZzKw1jtYI7q1ncGJYe14KJRvrmQLKtsHkwM6rsa2mms8htW1gLV3YVoW3XLGMsm3avn8PdU7X/JcyxrwM4GUAOH/+PHtW2qxQ0B5Ip2FWLuvPa2v6b6mkjaTansDBQe0RXVrSn1MpTV4CAa10lExqw+DECb3G6f31erWn8eJF/eB31uH4fDqS9HzjRw5RG+1WLGm1cP3q1Wb3pMdbJTuLi/X7Tg0Pa8n1rZKeYFAbt87eLI7x8c0L+R2N++BsJZnU+5ifrz9ujCYDllWd0ufxaAxZXNSkY6WmuXj9uiZ7L75Y3TjTKfTg9WrC4IyEbLW+pFkCNzam/x1mZupL5jtFCG7c2JxoNRLRv/Xv/V59VbNAQGPizZtbT8/bjYp7tHcOctvkUL8fs7Es4pn6jT8PD/Yh4N1583I+vjlx8rosBH1uFEo2MoUS7q5lYBsg5HdjZi2LYtng8cnItiuhubZYR7TVOQDIFsq4PB/Helpfb8DrwulDYQyFmsxdpn2l08nOLIAjNT8frhxrds2MiLgBRKCLAWkXOT2wjtu3tScW0J7oclkbTbatycnx49oIGB3Vxs36un7wu93aqLBt7Q0tFqvlrS1Lz0ejeuzkSe3xPXRIGw6Dg3odS6FSEz0RS0ql1iMQsVjz48kkcOFC/ftzZUWrg334w/p8d+7oaGokosnMjRvVNTXOSKrHo++1bBb47d+uVvM6dKg6JevIkeb30Ew6rddHo9UqiwMDem8rK5pg5CvT9n0+fb8vL9cnOo54XEc/rl6tXx9TKGg1uXBYv0ZG9Hc0Fhbo62s+6uH367qU2pE0QJOg/n6Nc/cr2wzo7z10SONZqaT34nT83C/ZoX2pJ+LJXrMqZaFn1jNYSeVhiWA80vfA+98UW0yLGwv7MT0awn95bxn9fR70+z0YCHggIlhK5LEYzm/7d45H/JiLbV7r43VbGNpiNMoYgzfurSOTr5auzBTKeHMmhvedGHqg5I52T6f/67wK4JSIHIcGjk8D+K6Ga74A4LsBfBXAdwD4E86J3X21i1eLxWqhgNlZnQvvzHFPpTTJsSxtLAQC2phaW6sWHjBGp3xMTOg1xaI2po4f1waDx6M9tuWyfk1O1vcoM9mhJnoilrhc+v9+s2pcrUZp7tzZPBKTyWingstV3/CfmwN+93d1XZ0zqpJO63vyIx/RtXm3b2snxfXrei4W04Tg8cd3Nk3UGWVxEhFHPq+/p/b15PPaWXLmTOvnu3mzeSEAZwraY4/pz88/r9PnnE6Y8XF93lZraKamdLra6qo+VzRaXYvUrDR0LdvW5MYZmaotkOLIZqvl8qlr9EQ82Q9cluDoUBBHhx7+Q3so6EMmn9l03O9xIeB1YTLaPEguJXPbTnYGgl5Mj4ZwcyW1EVe9bgtPHo5uWT1uNV2oS3Qctg3MxbKYHj2gZQ27REeTnco81x8C8AfQ8o6/Yox5R0R+GsAFY8wXAPwygH8vItcBrEGDDu0yZ1PPpSX9UDdGe2GvXNE3s1M5LRzWRkw4rD2iyWQ1mQH0Wo9HE6KVFf3+6FH9GhvTL0Abajdv6vm33tLEyJlywnU51KhXYomIroF7773N51qNDtQ2yI3R940z6nrtmr53T5/WzoZ4XJMfp3PBkcnotc5i+UhEp2Y5I7pnztRfvx0TE/qcjeto/H5NChpHqsbHqyM9zTgdHs76vpUVvbdoVDtJHB6Priva7toil6u+4mOtZpUoAf29776rSVappMmRs06nUSjERKfb9Eo86TVHhwJYTuaRK1aTChHgkbHQltPUBDsrcX1sOIjxqB/r6SIsCxgO+u5bJju/RcGFXPHBijHQ7un4uJsx5hUArzQc+8ma73MA/nqn74Pu75lnNPGYndUAMzenDalcrjr1plCoTkdzuXRNwLvvamPpyhVNVBIJbbDk89WvUqk6WjM3V23A9fdX96s4dar5vhGZjD5nIMBE6CDrlVjyyCPamHbWwnk8mui0SjYGBqrrYlZWqomOUxigWNRpa48/Xl3T46yPq1W7Vwug719nvdxWi/tbcbmA979f3//z8xoznLU1oZAmO7GYXjc0pO//SEQ7SBpHVDweHbn5i7/Qyoy1U93SaX38s88+eFJx5gzw6qubS2ufPt38+rfe0hFtRyql9+3xaEdP7dRdZ9So3dX0tsOZJuzx6H9LbuC+fb0ST3qJ3+PCc8cHdR1QtohsoQSPy8JSMo+RkA8et4Vik6RjrH/na2Z8bhcORba/SDHS13ph31bnaH/gJEPa4PFog+Lxx3X0ZmVFGzG5nDYMRDQRyWar3/f36/XBoD72yhXgT/9Uk6JgUBs/Xq82BE6f1sbdu+/qKM74uJboTSa1QeRsQJrP6xQcY3QzwNr9P4aHgfPnuZEodS8R4NFHNbnP5/V9tFVhgKkpTYwymfokYHS0OoUtm9XzTkdBs+eLRlvvF9Nsb6Dt8PuBp5/WL4dThGBgYPPIyciIrq15883qvff3azVIZy+ir3+9/jGBgL7fZ2aqe3/t1NgY8MIL9eW+T55sPrKTz29ODAFNJtxujWtvvKF/7/FxjY9f+Yomsa2Sp064eVPjrbM1QDCosZEdQtTNvG4Lx4eDeG8xiZVkHkAZQBEL8Rz8Xgu2JSjb1V6LiWgfRvsfbI3QToR8bhyK+LHQUEQh4HNh/AHXKNHuYbJDm3i92hB4+unqRqFer364O9PVnIaJ09i5eFEbTIcP6wduIlHdnG9oSJMU29bRI6cq0ltv6fQUZ4qKSHW6zblz2jCp7V0FtLH39tv6PETdzO2+/14tQHXvnuvXtYFbLuv7aWxM3yvOdDFnypffv3mUwe/XhCKR2DztTOTBk4hmRkf1/hoLETjHAa285qx1qV2fNz6uo1y1U9jGxjR5W1t7uPscGtre6Es2u3kTUkcgoAlNs5Gza9e0YEPj5p2dsLqqHUW10mkdvfrYxzjCQ90tlS/h7urmtTu5go1zkxEUyzZKtsFg0LuroyqPTfSj3+/BfDyLsjEYDvlwbCgId5M9gGh/YbJDTY2NaRLy4Q9r0rGwoGt44nGdtvHYY9pweeYZ7Uns69PpJ8ZowYHz5+uno5VK2pM7N6cNvFRq81QWp5GwuKjJTm2Z11rz89vb1I+oV/j9+p7w+bQ33+Gs/ykWNWkQAT7+8epIKaCdEufOadL03HNa3cxprPt8em6rks47JaJFBG7dqk5vGx/Xe63VrCCDz6dxxUmKGs9tVyKh8cbv1yRrJ43/YFD/fuXNa5HR319fwrqWMbrmsfF1dkKr2JjJaCLU7O9H1C1m1zIolu2mG4kmc0WcGtubYgAigqmhAKaGdqFHg9qKzUXasLamjaC+Pk12Tp/WUrDRaHUh89mzetypxuZwGijGAF/+cv0agLU1bficOaPTV1ZWmldBcub7O0lMsbj5GqBaIYnJDu132aw2ji1LOw+8O99nr87x49rx4IzmeL26UP/kSU1YIpFqpTHb1vdj7ZS2UAj40Id06mippNd3YnG9Zek9nTy5s8eNjmr8aZxuZ1nbG9UxRpO5ubnqsWBQk6/tVnj0eDRhuXGj/rjLpa+ncbS51m7FpFax8X7niPazRK6IK/NJ3FxOYWY9i6DPhSMDAfg91SB2v0ICRM2wuUiwbZ3+4GwOCmjD4H3v015hZ6fzQ4e23uEd0B5Up/c4FtMP3rt3tfHgNDaGh7X3cWBAk6hAQBMdZ98KZ78M53c3CoebV1Yi2k9u3ND1ac6UqLff1imfTlL/INxundI2N6fvIZ9Pp041a8hvlcS0cySnnSxLE5PXX9fRGUDf6866wPu5das+0QF0etfFi/p3266zZ/X33r6tHTODg9rJ09+v8enmzc2Pcbvrq8Z10shI8xEmy9qbQgl08GQLZWSLZQR9LvjcO9iNuIVCycbrd9ZRKhtEAx7MxbNI58u4uZzGmfEwLBEdJeb6GHoATHYI167VJzqANhDefFMX9bYqz9pKMAh88IPae3zzpvaINja8hoa09G0sVt0tHdCGoFOV6vRpnTpXuyeJZVX32yDarxIJ3Wumlm3rwvbh4Ycb4bEsbXA320SzF4TDOn02kdCpZNHo9qehNSssAOjociazs/U0J040LwceiejUv8uXq/sfeTw6pXe3RnampnSEqXHt0COPPPzoIdFWyrbB5bkEFhM6PcOytEjA6bHwluWh72chnkOprD1DbsvC0cEA7q5lUSjbiGeLGA758MihMDfvpAfC/2uo6egJoNPNnMpozdi2TqlxykJPTNR/2IfDW+9B4fFoo2ZlRaetRKP1U+OCQT1/+7YmRcGgjhDt115pIker95Rt6zqWo0d3937aZW1NO0bcbk22OjnC+iBVxZqts3E0bsz6MI4f13i3vKzxzSmisFtcLh2punNH78Hj0RG++428Ez2s9xaTG4kOoO+rmbUs/G4Xjg1vf2PRXLGM+XgOhZKNgYAHK6kcZmIZlMsGQZ8bg0Evzo67kcyXMDUYxNNT0aZreIi2g8kObdkIaHWuWNRyq85UE0AXTr///dU1A4A2Ahp7uB2HDmmP7VYf0H7/1ruuE+1HD/Ke2u/eeKN+5OTqVR3NGB/fu3tqNDbWfM+gYLA+LrWDU31yr7hcrUefiDrBts2m0suO2Vh228nOaiqPt2biGyWk35wp4M5KGgaAJYL1TBGr6QKmR4KI9nlxciTIRIceCv/vIYyNNT/uVFlr5r336hMdQEeB3nqr/lgopPuJNDpyhHPLqXe1ek+JtD63n83Pb54iZtu6Fmar0ZTdNj29eeTX5dI1P0T0cEq2qdvjplahvL1eHGMMLs8nNp7HGIPZda1IUqjZMDRbKGMlXcBA0IOh0ANuBEZUwZEdwiOP6FSy2h5Rt1urPLXSuAjYsbqqa2xq542fOaOjN3Nz2kAaH9eqS0S9anhY11U0lgg+fXp39mFpt/n55sdLJZ1GtVsL8+/H69X1gjMz1dLTU1Pd+Tcn2m+8bgtBnxvpfGnTuYHA9haLJbIl5IvVpCZdKKNsG4gIogEPwn4PElktKRjyuvHUkR0uGiZqgskOwefTcrSzs9pACAR05KWd8/G3u6EfUa948kld17GwoOs6JierFQepc1wuXRPVreuiiPazU2MhvHkvVrfxrsslODGyvSls0jCfqHZNr9uyMBntw2RUp5QMhbxwsdQ0tQGTHQKgDYSpqe3vUD4x0bz86tAQqwEROUZGemPR+Ph486ILHk9vvD4i2p7hkA/njw3i3loG2WIZIZ8bU4MBBH3ba072+z0IeF3IFHT+a8Djht/jQq5YxkDAU3fteKTFPHqiHWKyQw/kkUd0ylpt6VO/X6fvfPWrun5neFg34Wu17oeIusP4uI721m6oaVlaPv5BqpCtrek+ROm0rg08eVJLOhPR/hfp8yAy+eBv2McmI7h4L4ZiZY3OseEAYpkiBoPaUyoCHBkM4BD31KE2YbJDD8Tj0bnxi4vVjUGTSa3Q5Egmda7/Bz/ITUCJut1TT+nUsOVlXdM3MfFg7+uFBeDChepmq06ceOEF3byTiLqbMQa2QcspaJE+D16cHsZSUktPR/u86O9zYzVdQLFsYyDghd+zi7Xcqecx2aEHJqILkw8d0lLUjZXYAN19/NYt4NFHd//+iKi9BgZ2vslwoytXUDffH9DCJVevasJDRN3JGIMby2nMrGdQKhuE/W6cHA1huEk1NZclm6apNbuOqB1YepraIhZrvX/I2tru3gsR7U/Foo7kNMM4QbR/lW2D5WQeq6k87Bblp68uJnF7JY1SWc8ncyW8eS+GWKawm7dKtAlHdqgttprOwilsRATo+h63W0tWN2KcINqfFuI5XFlIbCQxXreFc5ORjTU2gO6RMxfLbnqsMcDdtQyi2yxNTdQJHNmhtgiHW5eWPnZsV2+FiPYpy2pdEppxgmj/yRRKeGcuvpHoAJrYvDkTQ6lmI1mdPesAAA5sSURBVNFssdxydkc6v492HqYDqSPJjoj8byJyRUTeEpHPi0jT3SVE5LaIvC0iF0XkQifuhXbPs8/Wl6H1enWvEe6vQw+D8aS3nDmjJe6d/TVcLq3GdvLk3t4XHQyMJzszH89tWmMHAOWywWIyv/FzwOuq2zOnVmibZamJOqVT/wd+CcCPGWNKIvJzAH4MwD9sce1HjTErHboP2gXz88D16zoXPxwGzp7VstPhMFoGP6IdYDzpArOzWk46ldJy0tPTWrykkWVpJ8ijjwKZDBAManVHol3CeLIDtSM6jco15zwuC4cHAri7mqm7xrKAqaFAx+6PaDs60hQ1xvyhMcaZlf01AIc78Xto783OahnZWAwol/Xfy5e1wcNEh9qB8WT/u3cPeP11LUNfLgPr68Crr2qZ6Va8XiAaZaJDu4vxZGdq1+U0GgjWv3lPjYZwcjQEn0c//KMBD546MoBIH9/ktLd2ozn6fQC+2OKcAfCHIvKaiLy01ZOIyEsickFELiwvL7f9JunBXLu2s+NED+mh4wljSfu9997OjhPtE4wn9zEc8mI4vLkk9ORAH8L++iRGRHB8OIgPnhrBJ86O4fyxwS2TJaLd8sDT2ETkjwA0maSAf2yM+Z3KNf8YQAnAf2jxNC8aY2ZFZBTAl0TkijHmz5pdaIx5GcDLAHD+/PnW46q0a4xpXUY2mdTz0nxPMaI6uxlPGEvaq1TS6WjNJBK7ey9EAONJO4kInpiMYCGRw3Iyr/vr9fsx2s/yidQ9HjjZMcZ8YqvzIvI9AL4NwMeNaba8DTDGzFb+XRKRzwN4DkDTZIf2HxGgrw/Ibq42ib4+Jjq0fYwn3cvlAnw+IJ/ffC4Y3P37IWI8aS/LEkxE+zAR7bv/xUT7UKeqsX0LgB8F8FeNMU37/EQkKCJh53sA3wTgUifuhzqnVQUlVlaidmE82d9EgBMnmp9jHKD9hvGE6ODp1JqdzwAIQ4d+L4rIZwFARCZE5JXKNWMA/kJE3gTwdQD/2Rjz+x26H+qQ48eBxx7Tnl1A/z17Vo8TtQnjyT43Pa3V1Zw44PcD585piWmifYbxZBcZY9Bi8Ixo13Sk9LQxZrrF8TkAn6x8fxPAk534/bS7TpzQ5KZY1MpKnL5G7cR40h2mp3Ukh3GA9jPGk91RKNm4tpTEYkL36RkK+XBqNIQg99yhPcDiwNQWIlpKlg0cooOLcYCIAODivRjmYznYthYrWknm8dqddRTL9l7fGh1ATHaIiIiIqC1WU3kkssVNxwslG/Ox3B7cER10THaIiIiIqC0yhXLLc6l8qeU5ok5hskNEREREbRHwulqeC3HNDu0BJjtERERE1BZDIR/6+zybjvs8Fsaj3IyUdh+THSIiIiJqm6enopiI9sFlCUSAkbAPzx4dgMfFZiftPo4nEhEREVHbeFwWzk704+xEP4wxEJZopD3EFJuIiIiIOoKJDu01JjtERERERNSTmOwQEREREVFPYrJDREREREQ9ickOERERERH1JCY7RERERETUk5jsEBERERFRT2KyQ0REREREPYnJDhERERER9SQmO0RERERE1JOY7BARERERUU/qWLIjIj8lIrMicrHy9ckW132LiFwVkesi8o86dT9E1J0YS4ioXRhPiA4ed4ef/18ZY/73VidFxAXgFwF8I4AZAK+KyBeMMZc7fF9E1F0YS4ioXRhPiA6QvZ7G9hyA68aYm8aYAoDfBPCpPb4nIuo+jCVE1C6MJ0Q9pNPJzg+JyFsi8isiMtDk/CSAezU/z1SObSIiL4nIBRG5sLy83Il7JaL9i7GEiNqF8YToAHmoZEdE/khELjX5+hSAfwvgJICnAMwD+PmH+V3GmJeNMeeNMedHRkYe5qmIaJ9hLCGidmE8IaJaD7Vmxxjzie1cJyL/F4Dfa3JqFsCRmp8PV44R0QHCWEJE7cJ4QkS1OlmNbbzmx/8GwKUml70K4JSIHBcRL4BPA/hCp+6JiLoPYwkRtQvjCdHB08lqbP9CRJ4CYADcBvB3AUBEJgD8kjHmk8aYkoj8EIA/AOAC8CvGmHc6eE9E1H0YS4ioXRhPiA6YjiU7xpi/1eL4HIBP1vz8CoBXOnUfRNTdGEuIqF0YT4gOnr0uPU1ERERERNQRTHaIiIiIiKgnMdkhIiIiIqKexGSHiIiIiIh6EpMdIiIiIiLqSUx2iIiIiIioJzHZISIiIiKinsRkh4iIiIiIehKTHSIiIiIi6klMdoiIiIiIqCcx2SEiIiIiop7EZIeIiIiIiHoSkx0iIiIiIupJTHaIiIiIiKgnMdkhIiIiIqKexGSHiIiIiIh6EpMdIiIiIiLqSUx2iIiIiIioJ7k78aQi8h8BnK78GAUQM8Y81eS62wCSAMoASsaY8524HyLqXownRNQujCdEB09Hkh1jzN9wvheRnwcQ3+LyjxpjVjpxH/9/e3cTKld9xnH8+yNBFyL0xRrfIdBQ0EWFhkDBQitiUje+UCWuLBXSha7cVMmigpu2tLgotZBSMZuahoIYGolGKYQu2jSCbRNN8KKW5tYa+rboxqI+Xdxzy5jeJDdz59yZ85/vBy73vM3M8+c/9wfPnDPnSho+80TSpJgn0vzppdlZliTAfcCtfb6OpPaZJ5ImxTyR5kff39n5EvBeVb15jv0FvJTk1SS7eq5F0rCZJ5ImxTyR5sTYZ3aSvAxctcKu3VX1fLd8P/DseZ7mlqpaTHIlcDjJyao6co7X2wXsArjhhhvGLVvSDFrPPDFLpLaZJ5JGpar6eeJkI7AIfKGqTq/i+MeBf1fV9y907NatW+vYsWNrL1LSipK8OktfyO0rT8wSqX/miaRJGSdP+ryM7Tbg5LmCJMllSS5fXgZuB473WI+k4TJPJE2KeSLNkT6bnZ2cdYo4yTVJXuhWNwG/TvJ74ChwsKoO9ViPpOEyTyRNinkizZHe7sZWVV9fYdtfgDu65beAz/f1+pLaYZ5ImhTzRJovfd+NTZIkSZKmwmZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1yWZHkiRJUpNsdiRJkiQ1aU3NTpJ7k5xI8lGSrWfteyzJQpJTSbaf4/Gbk/y2O+7nSS5ZSz2Shss8kTQp5omkZWs9s3McuAc4MroxyY3ATuAmYAfwVJINKzz+u8CTVfVZ4J/Ag2usR9JwmSeSJsU8kQSssdmpqjeq6tQKu+4E9lXV+1X1NrAAbBs9IEmAW4FfdJv2AnetpR5Jw2WeSJoU80TSso09Pe+1wG9G1k9320Z9GvhXVX1wnmP+J8kuYFe3+n6S4xOqdZZcAfxt2kX0oNVxQbtj+9y0Cxgx0TyZkyyBdt+bjmt4zJNha/m92erYWh0XjJEnF2x2krwMXLXCrt1V9fzFvuC4qmoPsKer6VhVbb3AQwbHcQ1Pq2NLcqyn5516nsxDlkC7Y3Ncw2OeDFur44J2x9bquGC8PLlgs1NVt41RyyJw/cj6dd22UX8HPpFkY/fpyUrHSGqIeSJpUswTSavR162nDwA7k1yaZDOwBTg6ekBVFfAr4GvdpgeAdTtTJGkwzBNJk2KeSHNmrbeevjvJaeCLwMEkLwJU1QlgP/A6cAh4qKo+7B7zQpJruqf4FvBIkgWWrpH96Spfes9a6p5hjmt4Wh3buo9rSnnS6vxBu2NzXMNjngxbq+OCdsfW6rhgjLFl6QMMSZIkSWpLX5exSZIkSdJU2exIkiRJatJgmp0k9yY5keSjJFvP2vdYkoUkp5Jsn1aNk5Dk8SSLSV7rfu6Ydk1rkWRHNy8LSR6ddj2TkuSdJH/s5qiX26qulyRPJzkz+v8hknwqyeEkb3a/PznNGifNPBkm82T2mSdt5olZMhyt5Mkks2QwzQ5wHLgHODK6McmNwE7gJmAH8FSSDetf3kQ9WVU3dz8vTLuYcXXz8CPgq8CNwP3dfLXiK90cDf1e9s+w9Lcz6lHglaraArzSrbfEPBkY82QwnsE8AZrME7NkOFrIk2eYUJYMptmpqjeq6tQKu+4E9lXV+1X1NrAAbFvf6nQO24CFqnqrqv4D7GNpvjRDquoI8I+zNt8J7O2W9wJ3rWtRPTNPBsk8GQDz5GPMk9lklgzAJLNkMM3OeVwL/Hlk/XS3bcgeTvKH7hTekE/3tzg3ywp4KcmrSXZNu5gebKqqd7vlvwKbplnMOmrxPWuezD7zpE2tvWfNkmFoOU/GypKN/dVz8ZK8DFy1wq7dVdXMP/Q63ziBHwNPsPRmfQL4AfCN9atOq3RLVS0muRI4nORk9ylEc6qqkgzuHvXmiXkyIObJjJuHPDFLmjEXeXIxWTJTzU5V3TbGwxaB60fWr+u2zazVjjPJT4Bf9lxOnwY3N6tVVYvd7zNJnmPptHhLYfJekqur6t0kVwNnpl3QxTJPPs48mV3myeybhzwxS9rQeJ6MlSUtXMZ2ANiZ5NIkm4EtwNEp1zS2bvKW3c3SFx+H6nfAliSbk1zC0hc1D0y5pjVLclmSy5eXgdsZ9jyt5ADwQLf8ANDEJ5erYJ7MLvNkuMyTgeeJWTIMc5AnY2XJTJ3ZOZ8kdwM/BD4DHEzyWlVtr6oTSfYDrwMfAA9V1YfTrHWNvpfkZpZOFb8DfHO65Yyvqj5I8jDwIrABeLqqTky5rEnYBDyXBJb+hn5WVYemW9L4kjwLfBm4Islp4NvAd4D9SR4E/gTcN70KJ888GR7zZBjMk2bzxCwZhmbyZJJZkqrBXTorSZIkSRfUwmVskiRJkvR/bHYkSZIkNclmR5IkSVKTbHYkSZIkNclmR5IkSVKTbHYkSZIkNclmR5IkSVKT/gtOkmn/neuLkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate many 2D (column) vectors \n",
    "np.random.seed(42)\n",
    "N = gaussian.rvs(0,1,(2,50))\n",
    "# copy and scale it\n",
    "S = np.copy(N)\n",
    "S[0,:] *= 4  # scale axis 0\n",
    "# rotate it\n",
    "f = +pi/4    # rotate by 45 degrees\n",
    "R = np.array([[cos(f), -sin(f)],\n",
    "           [sin(f),  cos(f)]]) \n",
    "X = R.dot(S)\n",
    "# shift it\n",
    "X += np.array([[1],[3]])\n",
    "\n",
    "plt.figure(figsize=(14,4)); #xlim(-10,10); ylim(-10,10);\n",
    "plt.subplot(1,3,1).set_aspect('equal'); plt.xlim(-10,10); plt.ylim(-10,10); plt.title('normal')\n",
    "plt.scatter(N[0,:],N[1,:], marker='o',color='b', s=50, alpha=0.3, edgecolor='none');\n",
    "plt.subplot(1,3,2).set_aspect('equal'); plt.xlim(-10,10); plt.ylim(-10,10); plt.title('scaled')\n",
    "plt.scatter(S[0,:],S[1,:], marker='o',color='b', s=50, alpha=0.3, edgecolor='none');\n",
    "plt.subplot(1,3,3).set_aspect('equal'); plt.xlim(-10,10); plt.ylim(-10,10); plt.title('final')\n",
    "plt.scatter(X[0,:],X[1,:],marker='o',s=50,alpha=0.3,edgecolor='none');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average\n",
      " [[0.34969052]\n",
      " [2.37483646]]\n",
      "Covariance\n",
      " [[6.99670167 6.59167706]\n",
      " [6.59167706 7.71554082]]\n"
     ]
    }
   ],
   "source": [
    "# subtract sample mean\n",
    "avg = np.mean(X, axis=1).reshape(X[:,1].size,1)\n",
    "X -= avg\n",
    "# sample covariance matrix\n",
    "C = X.dot(X.T) / (X[0,:].size-1) \n",
    "print (\"Average\\n\", avg)\n",
    "print (\"Covariance\\n\", C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.72610102, -0.68758803],\n",
       "        [ 0.68758803, -0.72610102]]), array([ 0.75465255, 13.95758994]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# eigen decomposition of sample covariance matrix\n",
    "L, E = np.linalg.eig(C)\n",
    "E, L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.68758803, -0.72610102],\n",
       "        [-0.72610102,  0.68758803]]), array([13.95758994,  0.75465255]))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# singular value decomposition of covariance yields the same\n",
    "E, L, E_same = np.linalg.svd(C)\n",
    "E, L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.0000000e+00, 2.7209678e-17],\n",
       "       [2.7209678e-17, 1.0000000e+00]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# it's a rotation!\n",
    "E.dot(E.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# also \n",
    "np.allclose( E.T, np.linalg.inv(E) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.68758803, -0.72610102],\n",
       "        [-0.72610102,  0.68758803]]), array([13.95758994,  0.75465255]))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# singular value decomposition of data matrix yields also the same\n",
    "U, W, V = np.linalg.svd(X)\n",
    "U, W**2 / (X[0,:].size-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.68758803, -0.72610102],\n",
       "        [-0.72610102,  0.68758803]]), array([13.95758994,  0.75465255]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# alternatively\n",
    "U, W**2 / (X.shape[1]-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[True, True]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check out the properties of U and V\n",
    "[ np.allclose( U.T.dot(U), np.eye(U[:,0].size) ), \n",
    "  np.allclose( V.T.dot(V), np.eye(V[:,0].size) )  ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# principle components from sklearn\n",
    "from sklearn import decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.68758803  0.72610102]\n",
      " [ 0.72610102 -0.68758803]] [13.95758994  0.75465255]\n"
     ]
    }
   ],
   "source": [
    "# object-oriented interface\n",
    "pca = decomposition.PCA(n_components=X[:,0].size)\n",
    "# sklearn uses a different convention\n",
    "pca.fit(X.T) # note the transpose\n",
    "#pca.transform(X.T)\n",
    "print (pca.components_.T, pca.explained_variance_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.68758803 -0.72610102]\n",
      " [-0.72610102  0.68758803]] [13.95758994  0.75465255]\n"
     ]
    }
   ],
   "source": [
    "# no more bug in the sklearn code - yeay! \n",
    "print (E, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.68758803 -0.72610102]\n",
      " [-0.72610102  0.68758803]] [13.67843814  0.7395595 ]\n"
     ]
    }
   ],
   "source": [
    "# use to yield results without bessel correction\n",
    "print (U, W**2 / X[0,:].size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAEzCAYAAACBoZBpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlspPd93/H3d27eN/fm7sq7kmzF8UXLdu0WTu04jpBGSeG0SoHWOYBt0hpo0BatUwOukaBA0zYtkDqNoCZG0yJN3CNOhESOLScp3ACNY8nRsdJKq5W0613uLu9rOMM5v/3jeUgNyRmez8PlLj8vgODM8zyc+fHhw888z+96zN0REZFoJO50AURE7iUKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCEUSqmb2JTObMLOLDcv6zexpM3st/N7X4mc/HW7zmpl9OoryiIjcKVGdqf4X4JPrln0W+GN3Pw/8cfh8DTPrB/4l8AHgYeBftgpfEZG7QSSh6u7fBGbWLX4U+M3w8W8CP9LkR38AeNrdZ9x9FniajeEsInLXiLNO9Yi73wof3waONNnmBHC94fmNcJmIyF0ptR9v4u5uZnsaD2tmF4ALAB0dHe978MEHIymbiMiKZ599dsrdh/byGnGG6riZHXP3W2Z2DJhoss0Y8NGG5yeB/9Psxdz9CeAJgNHRUX/mmWeiLa2IHHpmdm2vrxHn5f+TwEpr/qeB32+yzdeAT5hZX9hA9YlwmYjIXSmqLlW/Dfw/4AEzu2FmPw38a+D7zew14OPhc8xs1Mx+HcDdZ4BfBL4dfv1CuExE5K5kd+PUf7r8F5E4mNmz7j66l9fQiCoRkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQjFGqpm9oCZPdfwtWBmP7dum4+a2XzDNp+Ps0wiInFKxfni7v4q8G4AM0sCY8BXmmz6f939h+Isi4jIftjPy/+PAa+7+7V9fE8RkX21n6H6GPDbLdZ9yMyeN7OvmtlD+1gmEZFI7UuomlkG+GHgfzZZ/R3gtLu/C/iPwO+1eI0LZvaMmT0zOTkZX2FFRPZgv85UfxD4jruPr1/h7gvung8fPwWkzWywyXZPuPuou48ODQ3FX2IRkV3Yr1D9cVpc+pvZUTOz8PHDYZmm96lcIiKRirX1H8DMOoDvB/5+w7KfAXD3x4FPAT9rZlWgCDzm7h53uURE4hB7qLr7EjCwbtnjDY+/CHwx7nKIiOwHjagSEYmQQlVEJEIKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCClURUQipFAVEYmQQlVEJEIKVRGRCClURUQipFAVEYlQ7KFqZlfN7EUze87Mnmmy3szsV8zsipm9YGbvjbtMIiJxSe3T+3yfu0+1WPeDwPnw6wPAr4XfRUTuOgfh8v9R4L964M+BXjM7dqcLJSKyG/sRqg583cyeNbMLTdafAK43PL8RLhMRuevsx+X/R9x9zMyGgafN7BV3/+ZOXyQM5AsAIyMjUZdRRCQSsZ+puvtY+H0C+Arw8LpNxoBTDc9PhsvWv84T7j7q7qNDQ0NxFVdEZE9iDVUz6zCzrpXHwCeAi+s2exL4e2EvgA8C8+5+K85yiYjEJe7L/yPAV8xs5b3+u7v/kZn9DIC7Pw48BTwCXAEKwE/GXCYRkdjEGqru/gbwribLH2947MA/jLMcIiL75SB0qRIRuWcoVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQjFFqpmdsrM/tTMXjazl8zsHzXZ5qNmNm9mz4Vfn4+rPCIi+yEV42tXgX/i7t8xsy7gWTN72t1fXrfd/3X3H4qxHCIi+ya2M1V3v+Xu3wkfLwKXgBNxvZ+IyEGwL3WqZnYGeA/wrSarP2Rmz5vZV83sof0oj4hIXOK8/AfAzDqB/w38nLsvrFv9HeC0u+fN7BHg94DzLV7nAnABYGRkJMYSi4jsXqxnqmaWJgjU33L3312/3t0X3D0fPn4KSJvZYLPXcvcn3H3U3UeHhobiLLaIyK7F2fpvwG8Al9z937fY5mi4HWb2cFie6bjKJCIStzgv/z8M/F3gRTN7Llz2L4ARAHd/HPgU8LNmVgWKwGPu7jGWSUQkVrGFqrv/GWBbbPNF4ItxlUFEZL9pRJWISIQUqiIiEYq9S5WI3FmzS2XypSqd2RR9HZk7XZx7nkJV5B5VqdV5/vocc4XK6rKe9jTvPtVLOqmL1Lhoz8o9qVytM50vsVSq3umi3DFXJvJrAhVgvlDhtfH8HSrR4aAzVblrFMs1kgkjk9r8XOC18UWuzxao14Pn/Z0Z3nmi59Cdnd2eX266fHxhmbcf6yLsIi4RU6jKgTeVL3H59iKFcg0zGOzM8uCxLrKp5IZtx+aKXJsurFk2ky9z6dYC33uyd1fvXyhXWSrV6MymaMtsfM+DyN2pt+jy3Wq5REOhKgdavlTlhRtzq2ed7jC5WKJcq/P+M/0bth+bLTZ9ncnFEuVqfcuz3Ea1unNxbJ7JxdLqsqM9Od5xrJtE4mCf5ZkZA51ZphrKvqK/I6Oz1BgdrushueuMzRZXA7XRfKHCfLGyYXm52mRjgjCu1Jqva+XKRH5NoEJwSf3G1N1RJ3l+uHPDh0g6leD8ka47VKLDQWeqcqAVK7WW65YrNXra0muW9XWkuTW38Wdy6STtO7h0d3duzjc/6x2bW+bc8MEPpo5sig/eN8DNuSJL5SodmRTHe9t2dLYuO6dQlQOtO5dqegkL0JXbePieHexgKl+m0nDGagbnhjt3dMlbd6jVmtc9Vnd4xnsnZVIJzgx23OliHCr6yJID7URfG9n0xsP0WG+O9szGUG3PpPjA2X5G+tsASCSMB492cbQnt6P3TSaM3vZ003XqQC+bUahKrGp1Z75YYXmTy/jNZFNJ3n+mn+O9Qbh2ZFOcP9LJO451b/qek/ky7k61WufSrUVeujnPTidAOzfcSXJdg1QyabxtqHNXv8vdqlytc3VqiYtj87w5tUSpuru/5WGhy3/ZtkK5ysRCCTMY6so2PVNsdH2mwOuTearhZfRAZ4aHjveQSSWo153bC8uMzRXJl6pkkwlO9rUzMtC+4XVy6STvON46RNd78cYcb04uMZUvUa07bekk88UKPW1pTvZtfP1WetszPHy2n+uzBZZKNbpyKU72tW35e99LCuUqz16bpVR5q8rj2vQS7zvdR1eu+Zn8YXd4jg7Zk+9OF7g8vrj6/LXxPPcf6WoaghD0LX319uKaZdP5MhdvzjPYkeWlm/NcurUABr1h2F0eX6RYrnL+SNeuuyzlS1VeHV/bal+s1Lg6vUR/x85CFYLGngePbj/Q7zVXJvJrAhWgWnMuj+d53+m+O1Sqg02hKltaKlXXBOqKy+OLDHZlmp653WjRX/T1iTy3MkXG50vUHXCYWapQ9wIJg4s353hzeomhrhz3H+nc8dnQcqXGdH5jw5Y73JxrPsJIWptqsi8hmKSlVvcN1SOiOlXZhokWre8A4wvN15Va1KFO5UtUa76hq9SlWwvhOihV6swulXn22uyO62IzSSOZbP6P3qzBSzaXaNFjIpHYYgb6Q0xHmWxpN3e4adVCXnMnl06QbegrWarUWK7UqdeD7k8r/SirNefmXPMz3lbaMylG+ttJrDuy2zJJzqpr0Y4d62lruny4K3fgR5XdKbr8ly0Nd+d4Y3Kp6bqhrmzT5SP97dyeX94wwun8UBfJhDHYlWV+uYI7lGtOMhF0Y+pry6yZ+KRQ3tmZaiqZWJ08ZSZfplKv05FNMdiR5dzQwe+wf9C8baiDfKnK7FJ5dVl3W5r7tzEqq1iuMb1UIp1MMNiZPTRVBQrVGNTrTrlWJ5tK3BNjrDuzKe4b6tgQrPcNddCZbX4I5dJJHj7bz9XpJeYKFTKpBCd722jLJHnm6iyphJFJJrg1X6RcrdOTSzPUleV479ozo44Wr7+Ztw11kkokuDFbpFQNRl3dN9RJT4t+p/vJPehiVq07vW1pUgd85qxUMsH7TvcxX6iwWKrQkdneRNdXJha5Nl1g5SInnUrwrpM99Lbf+318FaoRcnfemFri+kyBas3JphOcGejgVP/OWpwPovuGOhnqyjIedqka7spu2YiUSyebtpyfG+7gG5cmABgZaKevPc1CscpAR3ZNHV4mleBEb/PLz82YGWcGO3Y1kqhaq5MvVcmkEpF3nVpYrvDijXmK4dl3Mmncf6RrV7/jfutpT2/7Q2lyscTVqbUzhVWqdV4cm+fDbxu856sNFKoRenNqiTcbzuZKlTqv3l4kmbANZ2B3o65cOpK+ibfmS5xa17Wpt63OcrVOKmnU3RnoyHKuyYQgcfrLa7P85Y05KtU6neEgg3ef6oukDPW68/z1uTXdk2o159LNBbpyKbrvoT6f4wvNe1mUKnVmC2UGOptXGd0rFKoRcXeut+hGdG26cE+E6opa3RmbLTJTKJMKPzD6tzl0s1qrs9BkdqlkIkF3LslHHxiOurjb8uzVGf701cnV54vLVWYLZQx4/9mBPb/+9FJ5Q3/PFTfninQfvXdCtVZv3bBZOwRzuSpUG5SqNWp139VlX7XuaybxgGCqucnFEsVKjbZMkhO9bS0bdnarXvd9vZyq1Z1nrs6wuPzWbUpuzy9z/kgnpwe2vtxOmJFMWNN/vFSLrlBxK1frPPvd2Q3Llyt1Lo/neefJXnLpYIarYrlG3X3Hdb2bTTtYbTFxy4pCucobk0tML5VJhx9ipwfaD2x9/WBXdsOUiRA2RKpO9XAoVWu8fHOB6XzQwtmeSXL/0S4Gd3CZkkoYbZnkan1ZpVbntfE85VqdzmySqcUSU4ulbYfPVlbm9SyUauTSSUb6mw/xjNrYbHFNoK54Y3KJ471tW96yJJEwjnTnmnaVulNn8/PFCpUWwZYvVSlV69TdeenmAvPhPZ/aM0kePNa97TP0YGJoaHaitlnDz3Klxrevzq5+YFcIRjktlas8dLxnW++9345157g9v7ymxwDA/Ue7DsUtbWL/Dc3sk2b2qpldMbPPNlmfNbMvh+u/ZWZn4i7Tes9fn18NVAi68bxwY478Dm4at9I4smJldnqzoEvSijemlvY8ddzEwjIXx+YplIIAX67UuDy+yLXp5t2etlIoV5nOl7Y1Ucb0UvPO/rW6M1soN1233v1HOhnoXBskR3tynI3gw2Y30kmjI9t8rtVsKkFHOslffnduNVAhOEaevz637cEJuXSS000+9Hra0xzrzlEoV3l9Ms/l8cU1I8JuzBY2XAFB8KFa3GF3s/2SSBjvHenlnSd7ONabY2SgnQ/c139XNMhFIdYzVTNLAr8KfD9wA/i2mT3p7i83bPbTwKy7nzOzx4BfAv52nOVqNF+oNK3jq9eDs7IHjm6/b+OJ3jZSCePadIE3pvJ0ZpMc6c6tadyp1ZyF5eq2z3CaubruHkwrrk0XGOnf/mVhtVbn4s2F1flKEwk40dvO/Udazz262ZlGuqHHfbla59Z8cXUikqM9udWfTSUTvGekj3ypSqEc3I9+fZVLoVxlarFMIhF0NI+zwaq3PcPZwQ7mCxvPWN93uo+ZYrlpgNXqzthccduzVp0b7qK7Lc3t+WWqdWewI8uJvjZuLyxz6dbC6lnsd6cLHO3J8dDxbhaaXBVAcMa7WKrs+z2zlis1JhdLuNNyiDIEJxlHunMc6d7ZlIv3grgv/x8Grrj7GwBm9jvAo0BjqD4KfCF8/L+AL5qZ+W6G8ezC8iZnZ7uZrm7lQMqmEk3rlYA9B0Sr2y6Xq3UqNSeT2l6ovnJ7cc0E0PV6MLNUWzrZsirhWE+u6V062zPJ1flH86VgZqPGM6xr0wXed7pvTQh0ZlNN+7lemchzdeqts+7XxvM8dKKb4a74/kHff6Yfw3h9Mritc1smyQfv6+edJ3u5PtP8Qwxo2fjUynBXbs3vUakFPUTWH+2355cZ7srSlm4dmputi8PNuSKv3F5Yvb3N5XF423CnRqqtE/fl/wngesPzG+Gyptu4exWYB/be3LpNzWaPX9HdtvsW2RN9zS91etvTLTvMb1erRpJMKkF6m4091VqdicXmXV/GNhkaOtAZdHVqHAbanknyvad6V89uX729uOGSdblS4/XJre/tNLNUXhOoEJwRvnRzIdYZ99szKc4MdDDcleXMQDtnBzowjEqtTk97GndnYjE4o3xxbD6ozy5X9zygYCacmKSZicUSp5oMuYWgHnY/p95brtTWBOqK1yfyLCxvvNI7zO6ahiozuwBcABgZGYnsddvD+/asbzjJpZN7qgMa7Mxy/5EuXp/Kr96Wo7c9zfec2HvjwpnBdl64Pr9x+UDH9i/96970hnoA5S3C68xgB8d725grlkknEvS2p1fft1b3DQ0UK1qduTdq1cexVnOm8uUdz+C/4XVazKw0sbDM5fFFsqnk6q2vJxdLvHxzgXed6qVUra+Z5WqhWKVaK/LRB/b2L7TVX6szm+JdJ3u5PJ5nqVRdrQ7ZSbVUo+szBcbmilRqdfrCao/t9GSYXCy1PF4mFpbvqX62exV3qI4BpxqenwyXNdvmhpmlgB5gev0LufsTwBMAo6OjkVYNvP1YF53ZFDfni1RrzkBncLDt9TJ9ZKCd4705FpeDETq7GXLZzHBXjneeDFrcl0pV2jJB6/9ORm6t3Aiv2dj6/m10e8mkEk0vx40glBeKFZIJo6ctvTpKajt5v9k96be6X/1m3ctuzhW5OrVEoVwjm05wur9jTRVHqz7GU/kSC8UKuXSS47055goVah4MMR3qynJjtsg79nBFM9CZJZW0pt2qVj5ABjqzfKgzy3KlRiphq0Nb63XHjG1/kF4eX+S7DfXxt+eXmcqX+MDZgS3rZjfb9Zt0Sz2U4g7VbwPnzewsQXg+Bvyddds8CXwa+H/Ap4A/2a/61BVmxshAPF2SUslELPc0Wqm7dfdd9Vd0d7pyKV6+tYARnEX3tGVIJY2zQ7uvI3ttIs/NuSJzYUt5KmmcHQjOhrZzljnUleVWk3lPEwladnEbmytyLQzM9kxy9Ux6xe35ZV6+ubD6vFSpr84Pu/I3X+n5UKvXWa7USacSZJIJ3Fm9Ffb6+lBgRz1EmkkmjIeO93BxbH61GmChWKEzl2K+WKE9k1xtDFrpK7syv+3MUjkcMpzj/iNdJCz4cJhcLJGw4Bg52deGmVGq1rgxu7FuuFpzrs8WtpwgZagry2sTG+t+IRiyLG+JNVTdvWpmnwG+BiSBL7n7S2b2C8Az7v4k8BvAfzOzK8AMQfDKNu02UF+4Mc/kYonhrixTi2UmFkt0ZFN8+NzwrluUJxaWuT5T4ERvG6VKnWKlRrXmXJ1Z4q/cN7CtVvLhrhxHe0obGsPOD3c1vXIYmytyqSEwC+XaaoCuBOvVFl3Nrk4vcao/CJ3etgxXp2YZXyithlswEUsHgx0ZLrfoY7qT2163MtSV5cPnBhlfKPLSzQUy6QQJM96cXOLq1BIPHe9Z/UCq1Oo8e212dfYv9+BDY3G5QjKRWNOTZa5QYb5Y4exgB9dmlpgvVulqcrXUrPfLem2ZJG8b6uTKxNp68VP97YdikpSdiL1O1d2fAp5at+zzDY+XgR+Luxzylql8ebV+sz2TYmTgrcNgL8MIb4VBmE4muP9IJ4vLQcf5XDrB24a33/H7e070cKwnx/RSmYQZR3tyLRv31jdqrS6fXloN1UK5dW+Jat1JJ422dJKJhkCF4Ay17tCWTXGke2OvBzM2zGGwW0EjY5JUIkHjZ4c7XLq9wGBnhlQywa25jdMpAlyfKZJO2prGK8f51pvTvDGZJ2nG6xN5cukk9w11kGn4W+S22YvgzGAHA50ZxheWcQ8+DBSoG901DVUSnZkWDUkA0/nSrnsnNNZ5mtma3hNb1YeuN9CZ3XLijXrdW3aAXxkYAdCRSTUdBZZNJ0iFdbDTSyXOD3cykS9RKNVIJY3BziwJC6oG3n6sm1TSuDW3HAxlziY5NxztdIKtemPUas5MocxwV65ldcNSuUoulaSxdmJqsczsUoX2dIqhrizduRQLy1W+O1PgXHjVsNMPhqgm1bmXKVQPoc3G2O9lfs/BzuyakWmN7xfHmO+VRqmx2aA1uz2bpL8jQyqRoL1hhNTZwQ5euLF5b4lytU42ndwQMO5QqTnZlPHg0W7uH+6iWvdYBiO0unUJgIX9BFp94KUTCXLrbhezMsJtpZvdyEA7Y7NF5osVyrU6PW1pzkf8wSAK1UPpWE+Oq1NLG+oIk0ljuCu7OoHHXKFCOhlM4LGdngUnetuYWCyt6VJlBg8c7Ypl1vfJxRIzhRKT4bDOuWKF6XyZc8OdnBl4ax7X4e6gt8SbU0vkl6u0Z5KcHuxY02Wupz3dtCdEJpWgveHyOJEwMjFNYDPcnW06sCKdSjAQNnYe681xbWZpw6CD0wPtVOr1Nd2eamHVxsoVQyqR4PRAB9V6ndHT/QyqgSkWCtVDqD2T4qHjPVy6vbDahzadCm5DUqv7mgk8litBZ/5ipbZlC/HKmO+JxRLT+TKZlHGspy2yrmSN3J3L44v0t2dxNyYWlimFZe7MpTZMzrLVkMmzgx1MLpY2dG26b6hj32YBG+7KcbK/zI2Zt7p3JRPG9xzvXi1DOpyJ/8pEnql8CSMYDnr+SCfzxQqv3FpcHQl4pDtHWzq54Qy4O5dWoMZIoXqI1OtB95nb88s4cKInR2dbikwySX97hkTCuDy+cTQUBBN7nB5oX+0Y38p+jfkulGur9akDHRkGOjKr3cuSu+gR0Z5J8fDZfq5NF5gvVsimEpzqb9/RTGVRePBoNyd625hZKpNKJhjuym5o4GvPpPjek70bfnawM8uHz2VYWK6SCG+g+OzV2TVn4IlEMFuUxEehGpFqrc5csULSbM0Io4PkhbH5NWP988tVutvSjJ7uWz0TatW9pl4Pts927u9481aaVSes7PPdVjW0Z1K8/djG27/st700BpkFAy5WvP9sP7fmloMPinRwe5o4rhxW1OvOxGIp/FAwjvXkDl3DlkI1AmNzRS6PL65eSrdlkrzzZM+BGro3VyivCdQVC8UKE4ul1X6QQfea5sG63a43+yGXTtLXkWk6JPZeusvCXqWTiX2ZZxeCQP3L63Nr/ibXZwo8eKz70Ez7B/swn+q9bmG5wqWbb9VNQjA7/PPX56gfoPF7c4XWHbznG85OT/W1Nx1O2t+ZieQMp1yt8+bUEhfH5rkykd/VTGArHjreTee6CXGO97ZxssVkNhKvm/PFDR9y7sHw2DgnwzlodKa6R81msIdgKOT0Ujny26fs1mZnmdmG7kE94aQvVybyFMs1zIJO3lFcFhfKVZ65Orum8/r12QLvHelbc8m6Xbl0kg/eN8DsUpnl8FbUUd8BVbZvqkl3Oljbz/Yw0BG4R5vdX2iz+xLtt6GuLJlUYsNonGTSONa79mA/0p1juCtLsVIjlUhE1ifz9YmlDe9fqzmvjS8yeqZ/168bx9wKsnObNRDupvHwbqXL/z1q9Q9txp5m949aMmG8Z6R3zfyx7Zkk7z7Z27RF38xoz6Qi7eQ+1eJWLHOFyoH6ANqt5UqN+ULlUF3qNmo1YU42nThQ/wtx05nqHh0Lb2I3v67O8vRA+4Fq2IGgVfkD9w1QKFepe+vROXFJmlFj45l9IrH5aKKDrlKr8/LNhdX5FJJJ48xAx6GbEX+oK8vpgXauNUwvuNL/+SD2homLQnWPgg7vfdycKzKZL5FKBBOAHOT6oztV73isJ7fmH27FcFculhFX++XSrYU1E3DXas7rE3na0sk9T6p9tzl/pIsTfUE/23QywWBn9q7+2+6GQjUCyYRxaoeTRB9G9w11ki9V18wP0NOe3nKk1kFWqtZa3tHgxmzh0IUqBB/ah7nB8PD+5rLvgnrdPhaWK6tj8O/2qePK1XrLWfGbTdEn9z6Fquy77lz6QA2M2IuOTIp0KtF0aO/d/oEhu6PWf5E9SCSMc8Mb72iQShpnBlUddBjpTFVkj070tpFNJbgxW2S5UqO3Pc3p/o5d35ZG7m4KVTk0ZpbKzBcr5NKJyHscDHZm931GKzmYFKpyz6vXneduzDHT0OvgSjrPe0b69r2vrtz7VKcq97xrM4U1gQrB3AyNt60WiYpCVe554wvNb6i3UKy0vHGgyG4pVOWet9mdXL3JsFmRvVCoyj2v1ZDhjuzhHvkj8VCoyj3vzED76h1FV6SSxjsOwK1T5N4Ty8e0mf1b4G8AZeB14Cfdfa7JdleBRaAGVN19NI7yyOGWSiYYPd3HZL60elO/Yz1tkU5rKLIirqPqaeB73P17gcvAz2+y7fe5+7sVqBKnRCK4y+v9R7o4PdCxo0B1d+YKZeYLFXyT+lkRiOlM1d2/3vD0z4FPxfE+InGbWSrz0s15SpVgbH8uneSh492624C0tB/XPz8FfLXFOge+bmbPmtmFfSiLyLaVqsENHFcCFYLZ/Z+7MXdP3KlA4rHrM1Uz+wZwtMmqz7n774fbfA6oAr/V4mU+4u5jZjYMPG1mr7j7N1u83wXgAsDIyMhuiy2ybePzJWpN7ohbqzm355c1f640tetQdfePb7bezH4C+CHgY96iIsrdx8LvE2b2FeBhoGmouvsTwBMAo6OjqtiS2JU3ORvdbJ0cbrFc/pvZJ4F/Bvywu2+8f0awTYeZda08Bj4BXIyjPCK70dfees7XPs2VKi3EVaf6RaCL4JL+OTN7HMDMjpvZU+E2R4A/M7Pngb8A/tDd/yim8ojs2EBnloHOjeE51JU9VHd+J1oOAAAKn0lEQVQHlZ2Jq/X/XIvlN4FHwsdvAO+K4/1FovKuk72MzRWZCO9DdaQ7y/GetjtcKjnINEZPZBMJ3dRRdkhDSkREIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQjFFqpm9gUzGzOz58KvR1ps90kze9XMrpjZZ+Mqj4jIfkjF/Pr/wd3/XauVZpYEfhX4fuAG8G0ze9LdX465XCIisbjTl/8PA1fc/Q13LwO/Azx6h8skIrJrcYfqZ8zsBTP7kpn1NVl/Arje8PxGuExE5K60p1A1s2+Y2cUmX48Cvwa8DXg3cAv45T2+1wUze8bMnpmcnNzLS4mIxGZPdaru/vHtbGdm/xn4gyarxoBTDc9PhsuavdcTwBMAo6OjvrOSiojsjzhb/481PP1R4GKTzb4NnDezs2aWAR4DnoyrTCIicYuz9f/fmNm7AQeuAn8fwMyOA7/u7o+4e9XMPgN8DUgCX3L3l2Isk4hIrGILVXf/uy2W3wQeaXj+FPBUXOUQEdlPd7pLlYjIPUWhKiISIYWqiEiEFKoiIhFSqIqIREihKiISIYWqiEiEFKoiIhFSqIqIREihKiISIYWqiEiEFKoiIhFSqIqIREihKiISIYWqiEiEFKoiIhFSqIqIREihKiISIYWqiEiEFKoiIhFSqIqIREihKiISIYWqiEiEFKoiIhFSqIqIRCgVx4ua2ZeBB8KnvcCcu7+7yXZXgUWgBlTdfTSO8oiI7JdYQtXd//bKYzP7ZWB+k82/z92n4iiHiMh+iyVUV5iZAX8L+Otxvo+IyEERd53qXwXG3f21Fusd+LqZPWtmF2Iui4hI7HZ9pmpm3wCONln1OXf//fDxjwO/vcnLfMTdx8xsGHjazF5x92+2eL8LwAWAkZGR3RZbRCRW5u7xvLBZChgD3ufuN7ax/ReAvLv/u622HR0d9WeeeWbvhRQRaWBmz+61wTzOy/+PA6+0ClQz6zCzrpXHwCeAizGWR0QkdnGG6mOsu/Q3s+Nm9lT49AjwZ2b2PPAXwB+6+x/FWB4RkdjF1vrv7j/RZNlN4JHw8RvAu+J6fxGRO0EjqkREIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIqRQFRGJkEJVRCRCClURkQgpVEVEIrSnUDWzHzOzl8ysbmaj69b9vJldMbNXzewHWvz8WTP7Vrjdl80ss5fyiIjcaXs9U70I/E3gm40LzewdwGPAQ8Angf9kZskmP/9LwH9w93PALPDTeyyPiMgdtadQdfdL7v5qk1WPAr/j7iV3fxO4AjzcuIGZGfDXgf8VLvpN4Ef2Uh4RkTstrjrVE8D1huc3wmWNBoA5d69uso2IyF0ltdUGZvYN4GiTVZ9z99+Pvkgty3EBuBA+LZnZxf167y0MAlN3uhAcnHKAytKKytLcQSrLA3t9gS1D1d0/vovXHQNONTw/GS5rNA30mlkqPFtttk1jOZ4AngAws2fcfbTVtvvpoJTloJQDVJZWVJbmDlpZ9voacV3+Pwk8ZmZZMzsLnAf+onEDd3fgT4FPhYs+Dezbma+ISBz22qXqR83sBvAh4A/N7GsA7v4S8D+Al4E/Av6hu9fCn3nKzI6HL/HPgX9sZlcI6lh/Yy/lERG507a8/N+Mu38F+EqLdf8K+FdNlj/S8PgN1vUK2KYndvEzcTkoZTko5QCVpRWVpbl7qiwWXIWLiEgUNExVRCRCBzZUD+IQ2PB1ngu/rprZcy22u2pmL4bb7bk1scV7fMHMxhrK80iL7T4Z7qcrZvbZmMryb83sFTN7wcy+Yma9LbaLbb9s9XuGjaZfDtd/y8zORPn+De9zysz+1MxeDo/ff9Rkm4+a2XzD3+7zcZQlfK9N97kFfiXcLy+Y2XtjKscDDb/vc2a2YGY/t26b2PaLmX3JzCYau2KaWb+ZPW1mr4Xf+1r87KfDbV4zs09v+WbufiC/gLcT9Bn7P8Bow/J3AM8DWeAs8DqQbPLz/wN4LHz8OPCzEZfvl4HPt1h3FRiMef98AfinW2yTDPfPfUAm3G/viKEsnwBS4eNfAn5pP/fLdn5P4B8Aj4ePHwO+HNPf5Rjw3vBxF3C5SVk+CvxBnMfHdvc58AjwVcCADwLf2ocyJYHbwOn92i/AXwPeC1xsWPZvgM+Gjz/b7LgF+oE3wu994eO+zd7rwJ6p+gEeAhu+/t8Cfjuq14zJw8AVd3/D3cvA7xDsv0i5+9f9rZFxf07Q53g/bef3fJTgOIDguPhY+HeMlLvfcvfvhI8XgUsc7JGCjwL/1QN/TtB3/FjM7/kx4HV3vxbz+6xy928CM+sWNx4TrTLiB4Cn3X3G3WeBpwnmM2npwIbqJg7CENi/Coy7+2st1jvwdTN7NhwJFpfPhJdsX2px6bKdfRW1nyI482kmrv2ynd9zdZvwuJgnOE5iE1YxvAf4VpPVHzKz583sq2b2UIzF2Gqf34lj5DFan5Ds134BOOLut8LHt4EjTbbZ8f7ZU5eqvbIDMgS20TbL9ONsfpb6EXcfM7Nh4GkzeyX8pIysLMCvAb9I8E/ziwTVET+10/eIoiwr+8XMPgdUgd9q8TKR7Je7gZl1Av8b+Dl3X1i3+jsEl775sC789wgGyMThQO3zsG3jh4Gfb7J6P/fLGu7uZhZJV6g7Gqp+QIbA7qRMZpYimO7wfZu8xlj4fcLMvkJwebrjA3m7+8fM/jPwB01WbWdfRVIWM/sJ4IeAj3lYGdXkNSLZL01s5/dc2eZG+DfsIThOImdmaYJA/S13/9316xtD1t2fMrP/ZGaD7h75+Pdt7PPIjpFt+kHgO+4+3qSs+7ZfQuNmdszdb4VVHhNNthkjqOtdcZKgnaelu/Hy/04Pgf048Iq732i20sw6zKxr5TFBI07kk7+sq/f60Rbv8W3gvAU9ITIEl11PxlCWTwL/DPhhdy+02CbO/bKd3/NJguMAguPiT1qF/16E9bS/AVxy93/fYpujK/W5ZvYwwf9h5AG/zX3+JPD3wl4AHwTmGy6J49DyKm+/9kuDxmOiVUZ8DfiEmfWFVWyfCJe1FkdLW0StdT9KUH9RAsaBrzWs+xxBa++rwA82LH8KOB4+vo8gbK8A/xPIRlSu/wL8zLplx4GnGt73+fDrJYLL4zj2z38DXgReCA+OY+vLEj5/hKAF+vUYy3KFoN7pufDr8fVliXu/NPs9gV8gCHqAXHgcXAmPi/ti2hcfIaiSeaFhfzwC/MzKcQN8JtwHzxM07P2VmMrSdJ+vK4sBvxrutxdp6GkTQ3k6CEKyp2HZvuwXgiC/BVTCXPlpgjr1PwZeA74B9IfbjgK/3vCzPxUeN1eAn9zqvTSiSkQkQnfj5b+IyIGlUBURiZBCVUQkQgpVEZEIKVRFRCKkUBURiZBCVUQkQgpVEZEI/X/8pBNAtBJDLAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# rotation\n",
    "A = E.T.dot(X);\n",
    "plt.figure(figsize=(5,5)); plt.xlim(-10,10); plt.ylim(-10,10);\n",
    "plt.scatter(A[0,:],A[1,:],marker='o',s=50,alpha=0.3,edgecolor='none');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVUAAAEzCAYAAACBoZBpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XlwnNd95vvvD92NlViJnSC4iItISiIlwZRkybJkybIk+1qTGSdRMjXjxJ7hJPe67vW9MzWTjO8kKadSlYwnmXKuM3EptivxlONlPJatsmVLsuOMLNtaSIqUuC/gBhALQew7Gn3uH6dBbN0gCLwv0ACeTxWK3W+/3X3QbD4871nNOYeIiAQja7kLICKymihURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCFEiomtlXzKzdzI5NOVZmZq+Y2dnkn6Vpnvvx5DlnzezjQZRHRGS5BFVT/VvgyRnHfg/4iXNuO/CT5P1pzKwM+EPgPmA/8IfpwldEZCUIJFSdc68CnTMOPwP8XfL23wH/JMVTPwS84pzrdM51Aa8wO5xFRFaMMNtUq5xzLcnbrUBVinM2AFem3G9KHhMRWZGiS/EmzjlnZouaD2tmB4ADAAUFBffefvvtgZRNRGTCoUOHOpxzFYt5jTBDtc3MapxzLWZWA7SnOKcZeGTK/TrgH1O9mHPuOeA5gIaGBnfw4MFgSysia56ZXVrsa4R5+f8CMNGb/3HgeynOeQl4wsxKkx1UTySPiYisSEENqfo68Etgp5k1mdkngT8FPmhmZ4HHk/cxswYz+xKAc64T+GPgreTPZ5PHRERWJFuJS//p8l9EwmBmh5xzDYt5Dc2oEhEJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJUKihamY7zezIlJ9eM/v0jHMeMbOeKef8QZhlEhEJUzTMF3fOnQb2AZhZBGgGnk9x6s+ccx8JsywiIkthKS//HwPOO+cuLeF7iogsqaUM1WeBr6d57AEzO2pmPzSzPUtYJhGRQC1JqJpZNvBR4H+kePgwsMk5txf4/4DvpnmNA2Z20MwOXrt2LbzCiogswlLVVJ8CDjvn2mY+4Jzrdc71J2+/CMTMrDzFec855xqccw0VFRXhl1hEZAGWKlR/gzSX/mZWbWaWvL0/WabrS1QuEZFAhdr7D2BmBcAHgX8z5djvADjnvgh8DPhdM4sDQ8CzzjkXdrlERMIQeqg65waA9TOOfXHK7S8AXwi7HCIiS0EzqkREAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJkEJVRCRAClURkQApVEVEAqRQFREJUOjTVEVk5YjH4dIlaG+HaBQ2boTq6uUu1cqiUBURAMbH4Re/gJ5uR05nC7ldLZwxo7+hlm0PKVnnS6EqIgBcvgw9PVDc+Da515tvHO98uZmR7Hpy9u9dxtKtHGpTFREAOjogu7djWqAC4KD/xGXo7l6egq0wClURASA7G7J7Um9VFI0C2sZoXhSqIgJAfT24yOwWwZwcKCoimaxyMwpVEQGgtBS2PVJHVnQyFnJzYfsOsEgWbNiwjKVbOfRfj4jcUL8zj9rfvpuBXxwl4uKsWwfEYnD33b59QG5KoSoi00TraymurfQ9VwAVFRCJLG+hVhCFqojMFo1q1P8CKVRFVpmxMWhqguFh305aVQV+E3hZCgpVkVWkqwveeMMH64SyMrj/fl3BLxWFqsgqcvTo9EAF6OyE8+dhR22/r8KOjfl2UlVhQ6FQFVkl+vuhry/1Y9ePNsGZI+CcP3DxIlRWwnveA1kaWRkkfZoiq0S6Smesp4Py174Lx4/7JaiGh/0D7e3Q3Jz6SbJgoddUzewi0AeMA3HnXMOMxw34PPA0MAj8lnPucNjlElltCgr8zKcrV+DqVRgagrKxNvZde5myaBv046uzHR2wezfk5UFrq1/fTwKzVDXVR51z+2YGatJTwPbkzwHgr5eoTCKrTm0tNDb6ZoDEyBiVJ/6R0avXiCZGJ08aH/dtq6BL/xBkQpvqM8BXnXMOeN3MSsysxjnXstwFE1lpWlvhrrtg8NI1Sk/9gg3jx8jNSzB4poliiqC83J/Y2+v/1NTTwC3Ff1MOeNnMDpnZgRSPbwCuTLnflDwmIreopwcijLO16xDryxw5eQaWxWBhNeNd3TAw4E+MRmHLFg3wD8FS1FQfcs41m1kl8IqZnXLOvXqrL5IM5AMA9fX1QZdRZFXIz4f4lXayxscgK8JY4XqyezugIJ+sTVshJ9u3Edx3H9xxx3IXd1UKvabqnGtO/tkOPA/sn3FKMzC1pbwueWzm6zznnGtwzjVUVFSEVVyRFW3rVjCXuHF/sHITo+tKKSlJrjRVUgIPPAB33rmMpVzdQq2pmlkBkOWc60vefgL47IzTXgA+ZWbfAO4DetSeKrIwmzfD+H0VdF/OYnwsQSQWoaRhO+vLh2FkBB580F/2S2jCvvyvAp73o6aIAn/vnPuRmf0OgHPui8CL+OFU5/BDqn475DKJrGq37com8S/uIH74HaLRiQ7+XJ+4mzff9Pkzd1Stq4OampALvYqYm5hhsYI0NDS4gwcPLncxRDJbb68f3B+P+9lTlZU3nZY6saPqzO2obrvND21d7czsUJqhn/OWCUOqRCQMRUXJfVDmr7k59f5+jY2+1SAvL6CyrWIa+SsiN0ysSz2Tc+kfk+kUqiJyw1w7pmg3lflRqIrIDRs3pm52zcvzTbJycwpVEbmhuBj27vV7/U0oKID9+7X06nypo0pEptm40U+66uz0Q6pKS5e7RCuLQlVEZolE/OYAcut0+S+SoQYHfY/7yMhyl0RuhWqqIhkmHoe33/bL+EWGB4iNDrDh9kJ236tBoiuBQlUkwxw7Bq3N4xQ3Hia3qxWAvtNG09UN1H14rxaWznD62xHJIPG4n9VUePn4jUAFMOfoercJzp5dxtLJfChURTJIPA6JeIK8602zHhsbAy5fXvpCyS1RqIpkkNxcWJc3jiXGZz1WWAiMjs5+kmQUhapIhtl1V4x4QfG0Y5EIbKhjco8pyVgKVZEMU10Nuz+2m7LyLNat89ND9+yB/MIo3H77chdPbkK9/yIZqGRbOSWfeBguXID+fr+E35Ytfs6oZDSFqkimKiz0+03LiqLLfxGRAKmmKrIcWlvhyhU/hqqiwu8dFdU/x9VAf4siS+30aThzZvJ+R4cf8f/ggwrWVUCX/yLzlEhAWxu0tPgK5oI0N8NLL/lQbWvzO+2B36TvypXAyirLR/8tiszDtWtw+PDk2PtoFO6802/fPF/u0GH6v/V93MXLxIryyasq8rXU22/3A1E7OnwPv6xoClWRmxgbg7femqxUgq+pHjniF3CezyinkWNnaPnL7zLeM0Du9T642kekqZeKfRuIdHRAVZU2gVoldPkvchMtLdMDdYJz0DR7iv5sY2M0P/8Wo8OO8ew8EhG/V8l4/xBdl3qhp8eft3FjcIWWZaNQFbmJsbGFPTZhvP063cncxIzh9RtIRH2tdPDaAB292RyzOznVXsbQ0OLLK8srtMt/M9sIfBWoAhzwnHPu8zPOeQT4HnAheeg7zrnPhlUmkYWYa1uROXcY7eiAc+fgUhM5HVdhPA6RKC6azVDlZtzwEC0DxYze9quMuGo4C42NcP/9UFYW+K8hSyTMNtU48G+dc4fNrBA4ZGavOOdOzDjvZ865j4RYDpFFKSryw0gvXpx+vLp6jlBta/MNsc4RycumMHuYEZeA+CguWUvtHsnn6h1PUFhafeNp4+N+keqHHw7lV5ElEFqoOudagJbk7T4zOwlsAGaGqkjGu/NOX2O9etUPraqp8TuOpnX6tG90BTCj/IHtNP/0DHEXZahyEyQSvFt7NyUP3jHrqT09MDwMQ0P+ZTo7ISfHB/vWrdoqOtMtSe+/mW0G7gbeSPHwA2Z2FLgK/Dvn3PGlKJPIraqu9j9TDQ9DV5dfB/XGVs59fdDe7pMwKb+8gM0f3UvnxV467nyYgk0VbLqSy/Dw7Pcxg4EBeOONyQ6ywUE4ccK/35494fx+EozQQ9XM1gH/E/i0c653xsOHgU3OuX4zexr4LrA9zescAA4A1NfXh1hiWWuGh33TZ0cHxGKwadP8xp8eO+abBCYqpOV0cG/sHbLHBnwC5uT4qmVuLgCxnCyq7qqi6gnfyz+a5WuiM1VW+gX+U404uHgRduzw5ZTMFGrvv5nF8IH6Nefcd2Y+7pzrdc71J2+/CMTMLOUqvM6555xzDc65hgptSC4BGRmBn/3Mr7DX1+cvtd9+G06enPt5ly/750wEatbIENGf/yMXDnb4A1VVfsm+06d9e8GEKYP7t22bXfMtKvILU/XOrH4kJRL+ZSVzhdn7b8CXgZPOub9Ic0410Oacc2a2Hx/y18Mqk8hMFy6Q8hK8sdFXMqdcwU8zdauoWH8XlW+8wLqWcwCMJfKJ7dgyOa+1p8c3yG7Z4pM0KSsL3vMe/3BPD+TlTY40KChIHaxmkJ+/0N9WlkKYl/8PAv8CeNfMjiSP/UegHsA590XgY8DvmlkcGAKedW7i/36R8HV1pT6eSPigS9e7PzE+1cZGKTn9BrHBnsnn9g/6Gurevb5Ha9s22LXLT0Wd8R4tLb7muW4drF8/+djWrX4hq5n/Gurq0ge9ZIYwe/9fA+bsp3TOfQH4QlhlELmZuQIq2RSaUkWFD8O8601kjY8xnrsOeq4RiSZnm8bjvi2hosLPlJoRqMPD8MtfTr+ULyiA977Xv29ZGTQ0wKlTvlkiGoX6ep/Nktk091/WpKtX/SX+1au+CaCmZkrvPT7UiorSP3/bNl+TzBr1bQejhWXkdLdSXTY0OeRpdNS/cHHxrOefODG7bXQg2b91zz3+/sRog7Exn8lZmv+4IihUZU25fBlefx3efdfXCGtr/SX+RC9+WZmvXN5999yvkzvay8N7xriSXUx8CLKzI1Tu2kVBb4tvU8jK8uk4kZAztLSkft2WFl+OqWNR1dO/sihUZc24dAneeQfOnvX3h4cnO6T27vWX2B/4gO8wSmtiAGlvL9mRCLdFo7B1LJl8USje6C/3Kypg//6l+LUkwyhUZUVraprclaSy0gdkuprduXO+c2hmb39LC5SX+8fmDNTBQfjSlxi7dJWheAzKSijcvQnLyfaX+QMDvopZW+sLMofqat/0kOq4ZkytbApVWbGOH/c1zQnd3T4gH3po9q4k4+M+E+Nx39RpNhm+Q0N++NLgoG8WqK2d3hMP+MD8y7+k99Uj9PWDi8QYKa6k/cIgtR+6i/xYDN7//nmXfc8eP7pgYGDyWH4+7N59a5+BZB6FqqxIQ0O+g2mmvj5fc525gH5Wlq8ZNjf7DqKODj+MqbraB+qpU34m1cWL/mdiFNQNP/wh/Y1t9CU7l2x8jJzOqwxFYpx4rZOC+AYG477CWlNz89pmbi488ogvU1+f3426tladUauBQlVWpM7O2WM4J1y/PjtUz571PegTnVHg+5Oamvyl/8aNfhLUhHPn/JjQwsLkC3Z00Duahx8l6N/YcIy099AyOERix3qyrvqQrK2Fe++9+e+QlXVr27HIyqBQlRVprvGlqR67dMmHpplvIigrm2y/vOOOWcNIAWhtilMYu+gn+be0wMgQo+tKye7vBHxTwvDQGP2VlYyX1DMxcOrqVT+mVLOp1yaFqqxI69f7wfJT2yTBh2Sq9XYmNuyrrPQ/iYSvKTY3p359Gxmm5OvPQW+jr9729VE4MMQIJQyX1hAd7GUoPk57+S7O7v4od5VM/6fU3q5QXasUqrIimcF998GhQ5NbPGVn+1pnirH2lJf7oJsw0Xa5dasP3KlNCTY2ys6//0+U9h+DWJbvQRoeJj8nlxwXpy9vPfF1JbRHN3Cw+jeo25ozq6arsaVrl0JVMtrAgJ9RVFQ0uxOnoMCvkN/b6y/FS0rSd/Tcfrtvh43HJ49NjM8fGPBX+IkEMD7Oppf/hrprR4gy7Fek6OuDsjIikSyqbi8ltvUu2sv3UFi/jdsv58wahmUGGzYE+SnISqJQlYw0NASHD/sgBN9OumdP6rCaazrphOJiH8CNjT6E163znVlFRb4WW13tF5SKnj1DpZ0kmmc+UMFXY3t6oKqKSF0tFb/6KBPLT27Y6ss5scBKJOKX7pvPttWyOilUJSO9+eb0pe9GRvw6pwUFvka6EAUFfluUVHKi49R3H4eDP4TWZj9odXR0stdrYoBrWZlP4aTKSvjgB+HaNZ+9FRWzx8jK2qJRcZJxOjtTryXqnO/FD8WhQ36w6tiYD86JRtGp7QUVFfDhD88ahBqJ+JpuTY0CVVRTlQw0MrKwxxYkkYC338Z9/wcM9idIdHWTS4RYVZVP90TCtxVUVMAf/mGKqVYi0ylUJeOUlvrKYKrB/aWls1dxWpSTJxl+9wwtjQnfLpooIrermfzy9ZRvX+/f6O674aMf1RgpmReFqmSc3Fw/1On8+dmPtbX5q3Tws6D27FnE8KW+Pjh0iAunHVlxw3CQlcXw+o0Mj4+TtWEjZfs2wRNP6Lpe5k3fFMlIu3f7nvnLl30zZ3Gxv33liv+zvx/eeguOHIF//s8XEKzvvgtnzzJy/CyRixAZ6iMRy8FFswFIRLNpy9lEWarVWUTmoG+LZKy6usm58WfO+M7306cnt252ztdmf/xjeOqpW3jhK1f8qinRKOPRHGCE8bxCbGyU8Zx8suKjDK2vY2j3g/Mbr5U0MuKHbF2/7gcNbNqUfo8rWb0UqrIi9PX5YUsTgTpVS3Kx/anboczpyhX/pxm5t20gcrGR8Ti4WDYD1VsZz11H1+0PsGVriqlZaUxsdT00NHmstdU3T9xkaVVZZTSkSlaEwsLpgTVVXt7sNQDmNCWZsyrKKbtvB/GCIhKxHEaLyuna9V5ya8tuKQzPn09dvtOnp4/KktVPNVVZETZt8u2qM7eULiryI55uepXe3u4HuY6N+QSOx2+0lZZuKSGnqoT2nhxy7n6UneVZbNx4a02pEzO/ZorH/WQsjcRaOxSqsiLk5MAzz8Df/71fYDorywfVxDqoc4bqa6/5a/OREV+trary6bx5843kzC8wNr/vDjbXLuziLTt77rLL2qFQlRWjvBw+8Qk4edK3r0ajviNr5845nnTuHLzySnK1FPywgYEB39BZUuLTODfXrxdYWLjgsm3a5Id7zbR+va9Jy9qhUJUVZd06eM975nmyc361k4lAnXq8rc2n8XyW6J+Hqio/DOzMmck21LKytDtUyyoWeqia2ZPA54EI8CXn3J/OeDwH+CpwL3Ad+HXn3MWwyyVLzznftNnZ6S+J6+rmvmxelGPH4Ac/8NtJj4/7doIpwwPcwCCjLptYIrh9oW67zddYe3r876ca6toUaqiaWQT4K+CDQBPwlpm94Jw7MeW0TwJdzrltZvYs8GfAr4dZLll6iYTPt46OyWOnT/uFpif2jArMiy/CV7/q3/T6dd9+eu2an25aVkZPD7T359LYWE+ixwfhrl3BTH2NRtUptdaFPaRqP3DOOdfonBsFvgE8M+OcZ4C/S97+NvCYmXY+X20uXJgeqOAvk48cCfiNBgbg29+eXCBg3Tr/098PFy/S15ug9VoWbTsfZrS4gnjcD4c6ceLmLy0yH2Ff/m8Arky53wTcl+4c51zczHqA9cCMf4KykrW0pD4+MOAvl1NtgXLLenvh61/305oKCiYD1Tk/jzUSoblkD8279zNYu+3G0/r74aWX/IiricWrN24MoDyyJq2YjiozOwAcAKhPtbObrFiBXJe88w585zt+tlRXl2+8LS726VhY6NPyttu4fP+vMxKf3FCqv9+PJnDOr0s9Pu5rz6Ojvo1U5FaFffnfDEz9P78ueSzlOWYWBYrxHVbTOOeec841OOcaKrQE24pTW5v6eEHBLU2vT210FF54wSdiRQU3No3q6YHubn/bDN7/forLpu/Q19LiAzUanb4oy7lzswcNiMxH2KH6FrDdzLaYWTbwLPDCjHNeAD6evP0x4B+cS7WSpqxkmzfPXlwkFvN9R/MxPOzbZZubU8z/P3VqcvXqWMxPuC8sJE6E4Wu9DMWjuEc/AE8+yfbt03v7J6a31tZOPz46mn5arMhcQr38T7aRfgp4CT+k6ivOueNm9lngoHPuBeDLwH83s3NAJz54ZZXJyvI9/R0dk0OqamtvvmRfIgGvvw6vvuoD0MxvXfLhD/te+5TPqa6lbbSUsaY2hkuqufSR/5PYlo3cN+RHGtx/P5w96yuxZWW+Yjvz4icS0UwoWRhbiZXChoYGd/DgweUuhiyBEyfg+ed9TXWqDRvg2WeTmwAOD8PnPuerl8D1Trie7OZseeCf0nubrw5XVPhAnaq11a/LOtPmzek3CZTVy8wOOecaFvMaWqVKMpZzfgz/zEA1l2Ds+Cmu/vi472nKzYWPfMRXL5ncNLC/7nZ6t+678bxr12bvcVVdDXv3+pcA/xJbt/oWBJGFWDG9/7L2TPTIT1XT9CYPvf45ivpbKfsZ8O0N8K/+FTz+uG8PePttul8fobNiB4M122a9Zqr1WOvr/SCB4WE/wysSmX2OyHwpVCVjZWX5qayNjT5gc4a6ePzV/0RsdJBIFLIMhs43w5/8BTlbt5G1dTM89hjRShi8OPv1CgshPz/1e5lNDhoQWQxd/ktG27dvcjjWnhPfIjY6iOFDcHgEBgdgsK2P03/+/RuX/du3zw7PSATuuGNJiy5rlGqqktHKy+E3PzbKia8fpeT6a6yjj0hBHgmLMnXOgF1v59AhePRR3z768MN+HkB3tw/Y+vr0tVSRIClUJSP19Phl9PouXWfTT/+WXZFuiusGobWLwf4+hvPLb+x86hxcLrqDdw/7IVp33OFHBWhvKFkOClXJOL298POfA/397Pjan7Cu6SSDkRiuvoKSoiKy2rqJ0sNYUQWJBLRlb+C1in+G656sne7bN7kTq8hSUqhKxjl3Dlz/AFtf+Dwl5w9hyfmi4ye7GN+5nfG8DlxrB2P5RVysvo+Xd38alxUhFvPTXp2D48dnz5ISWQoKVck43d1QdvI1YgNdOMsikUgwOpocDnX2OpH3PkC3lXL8if+bd8/kMJacabV58+TiLKOjvglh3ttWiwREoSoZJz8f8tov+i2jcwqJ93SRcGBAYixB38VuYk/ew979OVxPbjlVXj45gH/CzabAioRBF0eScbZuhfGcArAs2vM3MxL1A0ijMX85P1xWzfGNT1FdDY895ttOZwZqSYm2M5HloZqqZIbxcd/L1N5OZTSKe3Az/d85w0CkmP6quyhLXKc4NkRnxUbO/dr/i4vG6O/3awD09ExOEAA/yD+g/fxEbplCVZbfxFJUnZ03DlVlQ+WHtpP1vy4y0gdZVkdnSRXND/8mLhrDbLJ2unu3r91OrH6lPaJkOSlUZfk1NU0L1AlWVsq6//ARLvyknfHcAobLJ8dIVVdPv+TPzU2/EPZCdXT4NVyHhnxzwm23+dEFInNRqEpoEgm/qHRbm58mWlc3e91SwC8flYpzVOQPsO3DOzl1Chjxvfu1tXDXXWGW3Of8229P3u/pgatX4aGH1FYrc1OoSigSCXjzzel52dQEO3bAzp0zTp6rmz4Wo77KB/LgoF9FKjs7lCLf4Jzft2qmsTE/y+uee8J9f1nZ1PsvoWhpSV0BPXs2xTYl6bYuzc29UbXNyvI1xLADFXx4z1zDdUKKVgqRaRSqEor29tTHnUsRtqWlfpn9qQuZ5uXB/v3LMiUqFku/w6u2WJGb0eW/hCLdFX1uRxP57zTBlXGoqvLToKJR/+eGDXD9ur+/fn1Ae1ffuuxs327bPHPfX3wxReaiUJVQ1NX5nvOpCi+8Q1HXJT91tBN/Ld3SAg8+6GuksZjv1s8Ad93lh862tvr7kYjv/U/XUiEyQaEqoSgp8UvwnTjhO60iQ/0Ud19i2/YZ25V0d/sqYYalVTQK73mPb/8dGvITCjTtVeZDoSqh2bLF11ivXYPsluuU5aRpIr1+PeNCdUJenrZZkVujUF3lensnV78vL1/69x+73kv8eAvjnR2Mjg6RW5IioZaiS19kiShUV6lEAg4f9k2WE4qK4L77Zi8+EpYrPzlDy/86DQ5wCQYunKJ0VzVVe2smTzLze52IrBIaUrVKnT8/PVDB11qPHl2a9x88fZnub71Efksjsb7rgNG/YSedp68x2JUcBJqd7UfSL3KKknN+CNeVK34ZQJHlFEpN1cw+B/xvwChwHvht51x3ivMuAn3AOBB3zjWEUZ61qKkp9fH2dhgZCXm8ZWMj/d98iZxuP1g1u7eDsYLrDNRup2/zXbQX1bH5wU2+N2uR41AHBuCNN6aHaX29771fphFZssaFVVN9BbjDOXcXcAb4/TnOfdQ5t0+BGqzx8YU9dquc8533Bw/6n+bGEdyJ2XM8YwPdxPr8dKREfiGUlQUysP/w4dm108uXfa1VZDmEEqrOuZedc/Hk3dcBbcG2xCorUx8vLAx2q+a3355su2256jj1ymUaT49RuGn2+nuxAX+xUn5nzazHFqK/33fCpaJQleWyFG2qnwB+mOYxB7xsZofM7MASlGXN2LFjdnhGIn7saFA6OydnHeVcv0r5O/9A2clfEH/rMImrLRTvmh6eQ6NZDG3dQ+94QSC15bleIx5P/9jNTCx2LbIQC25TNbMfA6mmv3zGOfe95DmfAeLA19K8zEPOuWYzqwReMbNTzrlX07zfAeAAQL16i28qNxcefnhyy+b8fN/WGGQtdWIOf6yvk+LGw5hzjBUUQ1YWQ5c7qNpeSuFTd9J9sYumK9By91Nkldfz9tt+Faj77/c154UqKvK/Z6rFT6qqbv312trg9Gm/zF9urh9nu23bwssna9OCQ9U59/hcj5vZbwEfAR5zLvX//c655uSf7Wb2PLAfSBmqzrnngOcAGhoaVJeYh1jMr4gf5uvjHIUX3yEyMkgiOw+yIgzUbKcifha6u8mvr6erZgPyBW//AAARjUlEQVT9NTvJqp38z3B4GI4cgfe9b+Hvb+bXYTl0yA8hm1BY6KeU3oqODnjrrcla6vCwD/54HG6/feFllLUnrN7/J4F/D7zfOTeY5pwCIMs515e8/QTw2TDKI+HYEGun/d13KDpzkOjwAPHcAgartzJeUEThnftgoBt27OB86w4GRmYPju3u9svsLab2XF0N738/XLrkRzWUlflZXNFb/GafO5f6sv/CBdg+c2qtyBzCalP9AlCIv6Q/YmZfBDCzWjN7MXlOFfCamR0F3gR+4Jz7UUjlkaANDJDzzlvs3Dh0IxWjwwMUtZxh2zZHLDfi10LdtYvxWPrZBkG0X65bB3v2+CGvE4te3aq+vtTH4/H0a6uKpBJKTdU5l7Ilyjl3FXg6ebsR2BvG+8sSuHwZEglKS6HokRqG3uqExDh5eSNE6AZKfdtDLEZNjV+ceqaioszZ86mwMHV4RqNLNwNNVgfNqJKFmZJAkfxc1r1nF+vqSiEWo2swh7bKOxnbtgvw7ZvFxdOfHouFv8/Urdi2LfVkgS1bdOkvt0Zz/2X+envh1Cnfq9PV5a+Za2t9GuXn01WxncZeaC99P/H2IiKv+I6kjRv9hnmtrf5peXm+3TOT1lEpL4eGBt/739s72ft/qx1eIgpVmZ+BAfj5zycHgBYU+CaA0VHYsoWxMTh/DgZK64jnFwF+HOnRo77zqKDA52/Q20gHqbra/yQSy7KLi6wS+urI/Fy4MH1EfSTixxrFYpCTw/VEKd0b76R3675pT3Mu/ToEmUqBKouhmqrMra8PfvlLeO01f3/nTt/DBL4Xp64O9u1joK+SoROpX2Ixs5tEVhqFqqTX3g5f/rIfTNrW5qcanT8Pjz4KVVWMj/ta6NlEAT1xP3urrm72Cljp1iEQWY10oSPp/ehHftUS8Mv0mflq56FDgO/UuTRazVBWAdnZviVgYhbShNpaP1xVZK1QTXUV6e31NcexMR9kNTULXFN0aAjefRdeesknZH6+r27W1t7o+e/tHKM9bxt99XtuPK2+3rcMxGL+1JqajNkcVWTJKFRXiStXfE/7xAyly5f9MKH77rvFjpdEwrehDgz4zqh43F/+NzX56UqbNoFzXN//NH3nZo+KLynxQ6j2alqHrFG6/F8F4nE4dmz2lM+OjgX0vLe2Tq76vHkzACOj0NsZp+dKL/E4uF27GBjP5cIFuHhx9hTPTJklJbIcVFNdBTo60vewt7be4r56U5bRd/v20X62l3hHK+AYc2NcLqihreAZBhv9qYODvj+rrs5f8mdnax8/WdsUqqvAXNMob3mK5cRwKaCzN8albY8Rq7hGrK+T7tvv50LtQzQehX37/OiqS5f8LKnmZn//nnsya6aUyFJTqK4C69enX6x5w4ZbfLHKSj9Rv6eH635LKcaKKxiq3EzXzgfouuxrxb29vv102zbfDOucby2Ykskia5LaVFeBrCy4997ZNcStWxfQ+24GDzwAmzbhojFcJMpQ+Ua6dr0XIpGUNd+sLF8j1sIjIqqprhplZfD4436M/tiY7/lfcIdRcgmpdevv4szh6Q+tX+/bcGfWSCMRDZ8SAdVUV5VIxHcWbdoUTA98qgVQCgvhqaem10qjUV9TjsUW/54iK51qqmuNc3D1qq/SZmX5bvvy8pSnmvmw3LLF9/BnZ/uQnWi/bW/34VpVtbDV9kVWI/1TWEuc87vbtbVNHrtyxe9nvXNn2qeVlfmfqXJzNXRKJBWF6lrhHLzxBvz0p4yPxHHFxUTrN/jq59mzfhpUkPtXi6xRCtW14uhRRn72JtfODDI4CHCN3JM9VDyyh7zCGFy75htjRWRR1FG1FgwMkLh0hQuXI8lA9YZ7Rml84xrj46hRVCQg+pe0QnV3+6X3Ojv9+qWbN/txqelO7uqCgdz1FNI6/bH+ATp6YlRpPJRIIFRTXYH6+uAXv/C97/G4n4N//DicSLPyPnl5DA/DeG4Bg5WbcFPWA4znFtC9rUEj90UCoprqCnT+vN9Ub6YLF2D79hTjRcvKyK0qhuYeRkuqGF1XRmyoF5cVpXX/R9m7WXNLRYKiUF2BenpSH08k/EL9pfFrfs2/eNzP5d+4kbIn93O57Shjze0QjTFUXk/fpjvIry6ipmZpyy+ymoUWqmb2R8C/Bq4lD/1H59yLKc57Evg8EAG+5Jz707DKtFoUFPgFTWYyg/yLJ+D4Yd/xVFDg1/5rasIeeIA7PnkfZ94dobUpTjyngJoavyGqdg8VCU7YNdX/6pz7L+keNLMI8FfAB4Em4C0ze8E5l651UPAdUq2tsxel3jpykpzvfWuybSA/37cHdHZCczOxjRvZc08Oe+7Jmf2iIhKI5a6j7AfOOecanXOjwDeAZ5a5TBmvrAwaGmDdOn8/EoHbiq6xs/v16Y2tg4N+YD/4Xi0RCV3YNdVPmdm/BA4C/9Y51zXj8Q3AlSn3m4D7Qi7TqlBd7X9GR/2VftbblyE7xV/n4KAfHqDVTkSWxKJqqmb2YzM7luLnGeCvgduAfUAL8OeLfK8DZnbQzA5eu3bt5k9YI7Kzk22io6OTW5nONDbmp6GKSOgWVVN1zj0+n/PM7G+A76d4qBmY+q+9Lnks1Xs9BzwH0NDQ4FKds6ZNLHS6fbu/5B8b88cjEb+lamnp8pZPZI0Is/e/xjnXkrz7K8CxFKe9BWw3sy34MH0W+M2wyrSqbdkyuXXqvn1+eMD4uA/UHTuWt2wia0iYbar/2cz2AQ64CPwbADOrxQ+deto5FzezTwEv4YdUfcU5dzzEMq1esRi8731+BkBHh1/4tL7ej1MVkSVjbua4nBWgoaHBHTx4cLmLsSo451sKYjE/zlVkLTOzQ865hsW8hmZUrWGXL8OZMzA05Du8tmzxTbIKV5GFU6iuUc3NcPTo5P3RUb/qlXNzbgIgIjex3IP/ZZmcP5/6+IULfg0BEVkYheoaNXWx6qnGxnytVUQWRqG6RhUWpj6em+sXvRaRhVGorlE7dqTukNq2TR1VIouhUF2jKipg/36/OEsk4me47tvnRwCIyMKp938Nq6zU3ACRoKmmKiISIIWqiEiAFKoiIgFSqIqIBEihKiISIIWqiEiAFKoiIgFSqIqIBEihKiISIIWqiEiAFKoiIgFSqIqIBEihKiISIIWqiEiAFKoiIgFSqIqIBEihKiISoFBW/jezbwITu8eXAN3OuX0pzrsI9AHjQNw51xBGeURElkoooeqc+/WJ22b250DPHKc/6pzrCKMcIiJLLdQ9qszMgF8DPhDm+4iIZIqw21TfB7Q5586medwBL5vZITM7EHJZRERCt+Caqpn9GKhO8dBnnHPfS97+DeDrc7zMQ865ZjOrBF4xs1POuVfTvN8B4ABAfX39QostIhIqc86F88JmUaAZuNc51zSP8/8I6HfO/ZebndvQ0OAOHjy4+EKKiExhZocW22Ee5uX/48CpdIFqZgVmVjhxG3gCOBZieUREQhdmqD7LjEt/M6s1sxeTd6uA18zsKPAm8APn3I9CLI+ISOhC6/13zv1WimNXgaeTtxuBvWG9v4jIctCMKhGRAClURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAC0qVM3sV83suJklzKxhxmO/b2bnzOy0mX0ozfO3mNkbyfO+aWbZiymPiMhyW2xN9RjwT4FXpx40s93As8Ae4Engv5lZJMXz/wz4r865bUAX8MlFlkdEZFktKlSdcyedc6dTPPQM8A3n3Ihz7gJwDtg/9QQzM+ADwLeTh/4O+CeLKY+IyHILq011A3Blyv2m5LGp1gPdzrn4HOeIiKwo0ZudYGY/BqpTPPQZ59z3gi9S2nIcAA4k746Y2bGleu+bKAc6lrsQZE45QGVJR2VJLZPKsnOxL3DTUHXOPb6A120GNk65X5c8NtV1oMTMosnaaqpzppbjOeA5ADM76JxrSHfuUsqUsmRKOUBlSUdlSS3TyrLY1wjr8v8F4FkzyzGzLcB24M2pJzjnHPBT4GPJQx8HlqzmKyIShsUOqfoVM2sCHgB+YGYvATjnjgPfAk4APwL+D+fcePI5L5pZbfIl/gPw/5jZOXwb65cXUx4RkeV208v/uTjnngeeT/PYnwB/kuL401NuNzJjVMA8PbeA54QlU8qSKeUAlSUdlSW1VVUW81fhIiISBE1TFREJUMaGaiZOgU2+zpHkz0UzO5LmvItm9m7yvEX3JqZ5jz8ys+Yp5Xk6zXlPJj+nc2b2eyGV5XNmdsrM3jGz582sJM15oX0uN/s9k52m30w+/oaZbQ7y/ae8z0Yz+6mZnUh+f/+vFOc8YmY9U/7u/iCMsiTfa87P3Ly/TH4u75jZPSGVY+eU3/eImfWa2adnnBPa52JmXzGz9qlDMc2szMxeMbOzyT9L0zz348lzzprZx2/6Zs65jPwBduHHjP0j0DDl+G7gKJADbAHOA5EUz/8W8Gzy9heB3w24fH8O/EGaxy4C5SF/Pn8E/LubnBNJfj5bgezk57Y7hLI8AUSTt/8M+LOl/Fzm83sC/zvwxeTtZ4FvhvT3UgPck7xdCJxJUZZHgO+H+f2Y72cOPA38EDDgfuCNJShTBGgFNi3V5wI8DNwDHJty7D8Dv5e8/XupvrdAGdCY/LM0ebt0rvfK2Jqqy+ApsMnX/zXg60G9Zkj2A+ecc43OuVHgG/jPL1DOuZfd5My41/FjjpfSfH7PZ/DfA/Dfi8eSf4+Bcs61OOcOJ2/3ASfJ7JmCzwBfdd7r+LHjNSG/52PAeefcpZDf5wbn3KtA54zDU78T6TLiQ8ArzrlO51wX8Ap+PZO0MjZU55AJU2DfB7Q5586medwBL5vZoeRMsLB8KnnJ9pU0ly7z+ayC9gl8zSeVsD6X+fyeN85Jfi968N+T0CSbGO4G3kjx8ANmdtTMfmhme0Isxs0+8+X4jjxL+grJUn0uAFXOuZbk7VagKsU5t/z5LGpI1WJZhkyBnWqeZfoN5q6lPuScazazSuAVMzuV/J8ysLIAfw38Mf4fzR/jmyM+cavvEURZJj4XM/sMEAe+luZlAvlcVgIzWwf8T+DTzrneGQ8fxl/69ifbwr+LnyAThoz6zJN9Gx8Ffj/Fw0v5uUzjnHNmFshQqGUNVZchU2BvpUxmFsUvd3jvHK/RnPyz3cyex1+e3vIXeb6fj5n9DfD9FA/N57MKpCxm9lvAR4DHXLIxKsVrBPK5pDCf33PinKbk32Ex/nsSODOL4QP1a86578x8fGrIOudeNLP/ZmblzrnA57/P4zMP7DsyT08Bh51zbSnKumSfS1KbmdU451qSTR7tKc5pxrf1TqjD9/OktRIv/5d7CuzjwCnnXFOqB82swMwKJ27jO3ECX/xlRrvXr6R5j7eA7eZHQmTjL7teCKEsTwL/Hvioc24wzTlhfi7z+T1fwH8PwH8v/iFd+C9Gsp32y8BJ59xfpDmneqI918z24/8dBh7w8/zMXwD+ZXIUwP1Az5RL4jCkvcpbqs9liqnfiXQZ8RLwhJmVJpvYnkgeSy+MnraAeut+Bd9+MQK0AS9Neewz+N7e08BTU46/CNQmb2/Fh+054H8AOQGV62+B35lxrBZ4ccr7Hk3+HMdfHofx+fx34F3gneSXo2ZmWZL3n8b3QJ8PsSzn8O1OR5I/X5xZlrA/l1S/J/BZfNAD5Ca/B+eS34utIX0WD+GbZN6Z8nk8DfzOxPcG+FTyMziK79h7b0hlSfmZzyiLAX+V/NzeZcpImxDKU4APyeIpx5bkc8EHeQswlsyVT+Lb1H8CnAV+DJQlz20AvjTluZ9Ifm/OAb99s/fSjCoRkQCtxMt/EZGMpVAVEQmQQlVEJEAKVRGRAClURUQCpFAVEQmQQlVEJEAKVRGRAP3/SCxsNbxC8A0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# projected points\n",
    "P = E[:,0:1].dot(E[:,0:1].T).dot(X)\n",
    "R = X - P  # residuals\n",
    "\n",
    "plt.figure(figsize=(5,5)); plt.xlim(-10,10); plt.ylim(-10,10);\n",
    "# plt.quiver(P[0,:],P[1,:],R[0,:],R[1,:],angles='xy',scale_units='xy',scale=1,alpha=.5)\n",
    "plt.scatter(X[0,:],X[1,:],marker='o',color='b', s=50, alpha=0.3, edgecolor='none');\n",
    "plt.scatter(P[0,:],P[1,:],marker='o',color='r', s=50, alpha=0.3, edgecolor='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Whitening\n",
    "- Frequently used to preprocess data, e.g., in signal processing\n",
    "\n",
    ">$ \\displaystyle Z = \\Lambda^{-1/2}\\ E^T\\ X$\n",
    "\n",
    "- There are two things we are trying to accomplish with whitening:\n",
    "\n",
    "> Make the features less correlated with one another.\n",
    "\n",
    "> Give all of the features the same variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT4AAAEyCAYAAABj+rxLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHk5JREFUeJzt3etvnGeZBvDrnvH5bNeOnfjc5tA2h7bUaktb6EILdKGCrywCCfGhXxYJJFZoof/ASkjASiChil1pJSqhlQCxQpzKLv1AW0qTtGmasxM7B+dkx4njjOPj3Pvhzsuc3hmPZ96Zd2ae6ydFrWfGM4/T5spzvB9RVRARuSQSdgOIiMqNwUdEzmHwEZFzGHxE5BwGHxE5h8FHRM5h8BGRcxh8ROQcBh8ROacujA/t7e3VsbGxMD6aiGrYoUOH5lS1b7PXhRJ8Y2NjOHjwYBgfTUQ1TETO5/M6DnWJyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icw+AjIucw+IjIOQw+InIOg4+InMPgIyLnMPiIyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icw+AjIucw+IjIOQw+InJOYMEnIlEReU9EfhPUexIRlUKQPb5vADgR4PsREZVEIMEnIkMAPgfgp0G8HxFRKQXV4/shgG8DiGd7gYi8LCIHReTg7OxsQB9LRLR1RQefiLwE4LqqHsr1OlV9VVUnVHWir6+v2I8lIipYED2+ZwB8XkSmAfwcwCdF5GcBvC8RUUkUHXyq+h1VHVLVMQBfBPB/qvrloltGRFQi3MdHRM6pC/LNVPUNAG8E+Z5EREFjj4+InMPgIyLnMPiIyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icw+AjIucw+IjIOQw+InIOg4+InMPgIyLnMPiIyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icw+AjIucw+IjIOQw+InIOg4+InMPgIyLnMPiIyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icw+AjIucw+IjIOQw+InIOg4+InMPgIyLnMPiIyDkMPiJyDoOPiJzD4CMi5zD4iMg5DD4icg6Dj4icU3TwiciwiPxZRI6LyDER+UYQDSMiKpW6AN5jHcC3VPWwiLQDOCQir6vq8QDem6hirK0BIkBdEH9qKFRF/ydU1SsArtz790UROQFgEACDj2rC4iJw9Chw44YF37ZtwP79QHNz2C2jQgU6xyciYwAeA/BOkO9LFJa1NeCttyz0AEAVuHYNePttIB4Pt21UuMCCT0TaAPwCwDdV9bbP8y+LyEEROTg7OxvUxxKV1MWLwOpq5uOxmAUgVadAgk9E6mGh95qq/tLvNar6qqpOqOpEX19fEB9LVHJLS9mfi8XK1w4KVhCrugLgPwCcUNXvF98kIrO8DNy5E24bOjqyP9fZWb52ULCCWJ96BsBXABwVkffvPfZdVf1tAO9NDrp7F3j/fWBuzr5ubQX27bNFhXIbHAQmJzN7d93dAAcu1SuIVd2/AJAA2kIEAHjnHVtJ9cRiwLvvAs89B7S1lbct0Sjw9NPAqVPAlStAJGJhuHt3edtBweKOJArV7Cxw/TpQX2+Bcvduauh54nHg/Hlg797yt7GpCXjkEftFtYHBR6FQBQ4dsl6U5/RpoL8/+/csL5e+XeQGBh+F4vLl1NADLAwvXrThZDSa+T1dXZmPLS4mtpz09gI7dtj3E+XC4KNQpIeep77ehpbpvbuWFmBkJPWxmRngvfcsMAELwPPngaee8g9OIg+Dj0KRq1e2f7/tn7t0CVhft9XcXbssFD0bG3aMzAs9z/y8hd/995em3VQbGHwUih07rMeWrr7egi4SyR1e8/N2nMzPtWsMPsqNsyEUioEBYHQ09bFoFPjIR/Kbo8s1lOUwlzbDHh+F5sABYGwssZ1lx47U4Wwu3d027+d3pGxoKNBmUg1i8FGoOjpyHwvLRgSYmLDNzisricfHxy1AiXJh8FHV6uwEXnjB5vS87SytrWG3iqoBg4+qWiQCbN8ediuo2nBxg4icw+AjIudwqEsV6/Ztq8fX1lbYAghRNgw+qjgbG1bAILm0+7ZttorLPXoUBA51qeKcOJF5n8X16/Z4JVtYAI4fB44dS1xORJWJPT4qiYUFO2u7sWG9tf5+23unav/M5eJF++famu3Ra2y0jc0XL1ol5ko0OZkazOfO2cmUAwfCaxNlx+CjwE1NAR9+mPj6/HkLwNZWq7rS0QHs2WPH1tKpWuBduGC9PC8oe3stSPIJznJbWgJOnsx8/Px5K656333lbxPlxuCjgp07B0xPW5j19FiYtbTYcC/Z5csWhqOjtufu9m0rJf/kk5n3aIhYLy95qKtqlZr7+iov9ADg6tXMKjHJzzH4Kg+Djwpy8iRw5kzi69lZm9caH0+9aHtlxRYqFhbs8qChIfvV32/DQ78LhOrrbRFjYyPxWDQKNDRUZo8vV1EFFkWtTAw+2rL1devtpYvHE/NznjNnLPQAC4GNDRsCNjZakPmpq7O5vOvX7Q6O5mYLyIYG+4xKW9ndvt0WNJID38Nzw5WJfx/RlsViqb2xZNFoosLKyoq91guq9nY7Uzs7a5WTl5b8w6K724JxeNhuMxsetq87Oysv9ABr26OPpvbuRICHHuLdu5WKPT7asuZm+0PuF1odHVZq6tAh23wsYkEWj1tYzszYcHVpye7LePvtzFLxDz4I/PWvqe8vYo9XqsFBW4C5etXaPTBgv09UmRh8tGUNDTZPd+FC6uMiwAMP2GT+pz5lIbeyYtcyzs0Bb75podfQYD257m7/UvH33Qc8+yxw9qyFY1ubPd/dXd6fc6saGzOLq1JlYvBRQfbvt7m4Cxdszq+tzXpk3gpmXZ2FwMc/bvNf7e3WE1S1oXBy782vVHxnp1VjJioFBh8VJBKxy70fesiCL9tCxf332xaXDz6wkGxrs8WAxsbU9yIqJwYfFSUSyR56noEB275SV2eLHelYKp7KjX/XUlmIAI8/ntrTA2w4PDgYTpvIXezxUdkkl4pfWbFV0La2sFtFLmLwUVm5XCp+YcHmQ7u6KnM/oksYfEQlFosBBw/aGWXAVrX37rWN2RQOBh+V3Pq6HWW7ccPm+EZG3DnRoGpXYCYv6qytAUeO2Bafrq7w2uYyBl8NuX3bTkd0dlbOFpG1NeCttxK9HcA2LD/6aHlXc9fXbTvNjRu2Cj0yUp4N0Tdu+K9kq1p7GHzhYPDVgMVF4PDhRLg0NdkGY796d+U2PZ0aeoD9oT92zA7wlyOg19ft1EhyOy5csCKhpT5psbqa/bnki9CpvCqkX0CFisdtKJX8h3p5OXFWNmzXr/s/vroK3LxZ+PuurdkxOL/eVLqpqczwBaxu4Pp64W3IR09P9jJarNMXHvb4qpxXuimdVyLqoYfK36ZkdTn+D8v1XC6nTwOnTlmPqb7eNkd/5CPZN1LPzvo/vr5u4dvXV1g78tHUZOeXJydTH29vt+E2hYPBV+VyDZeWl8vXjmyGhvx7fR0dhS1wXLoEvP66hdnGhgXf9u02ZH7iCf/v8cpkbfW5oHjlqS5etJ7qtm1WsLXQ4Kfi8be+yuUaLlXCUGpwELh1y4abXnn21lY7xVGIN96w0k8e736Oujqbs2tqyvye4eHU7/GUc1V1xw4WJa0kDL4q19ZmQ6b0ElGdneU5Cra+bqu0s7MWPsPDNvRMtnevFSuYn7fhaG9vYeXj43Hr8fm5etXmDf2Cb2DAymBNTiZq/LW22j295CYGXw04cMC2ZnjXOfb3W9CU+nTAxoZtVfFKywPAlSt26dDu3amvbW4uPohXV63Si9+CRjye+/jbnj1WFsvbS5hr0YFqH4OvBohYr29kxObTrl2zBYChIZtLK5ULF1JDz3PmjG0TSS9IUKzGRpsbu3XLhrjJ9u7dfGtMY6P/cHN93X6/ZmYS1ZP37PHvPVJtYPBVgbU16+14Jd/9qNpevsuXE4+dPZsYZpbC3Jz/4/G49ayCntMSsd6td/3knTu2ODEwAHz0oxaITU1bD6y//c3a6/E2Oj/3HM/U1ioGXwXb2ACOHk30RBobrSfit+n22rXU0POcOGFDTL/eVzxuiw4zM/ZZAwPAzp35r3Tmet1mNfoKNTpqn3vunN3b0dFh4XT4sPXcRGyV95FH8ls1nZtLDT1PLGZTBywlX5sYfBXsgw9SJ/NXVuyxxsbMUxnJF3Ani8dt+Ot3IP6991LDcnLS3udjH8uvpzMyknmdJGDzcKVcUU5eIb10yX4Oj6r9TJEI8Nhjm7+X38Zmj98wnmoDT25UqNVV64n5mZrKfGyrl1ovLPj3EBcXs39uup4eOxqXHJKtrbafrlwLB9PT/o9fvpw5D+inpSX7c62tBTWJqgB7fBVqeTmx7y2d30mNwUH/EKiry9xeAuQ+LjY/n/+pgrExW0SZn7fP6unJ7/uCkm2Tdjxuf3lkG47HYrYCDVi704+u1dezbFQtCyT4RORFAP8OIArgp6r6b0G8r8taW+0Pn1+vxW/TbU+Pzf+dPp0IzGjUjnL5zXXlWgBIvw82FrOFkps37bnx8dRjXnV1dhohDD09/j3Upqbs99pOTtrcp2d11X6v6urs987ryZZqnpLCV3TwiUgUwI8BfArAJQDvisj/qOrxYt/bZdEosGuXHaRPf3znTv/v2b07cUQsGrV5wGw9nv5+G+YtLWW+f3JP584d4C9/SQTw7ds2D/jII7l7hbGY/WpvL+3F2rt2WXvSe2x79vgP8RcXU0MPsIATsbnN5mYGnguC6PE9AWBSVc8BgIj8HMAXADD4ivTAA9ZzmZqyIV13t/1Bz7U3r6XFhp+bEQGeesoWBrxhb2ur9XSS573OnPHvdZ48aSGbHi4bG/ae3jBSxF534EDitQsLVmTA20w8Nlb4lpv2dgus9B5pth6o37wmYD29uTn7PS/GhQu24hyL2emZ3bvD6w1TdkEE3yCA5LW9SwCeTH+RiLwM4GUAGGFZirwNDpbu6FlrK/Dss9br29iwEEk3P+//vSsriR5dsuPHE6EHWKBcvGiBtGeP9bjefNM+D7Ce2rFjFuwPP1zYz9HWZj3QfGSbN93suXycO2c/i+fmTdsj+OSTpa0AQ1tXtlVdVX1VVSdUdaKP/xdUlJYW/9ADss8FimTuDfRKYfnxzhKfPZsIvWRTU/mtwiZbW0ucvc1Xrk3VxVyCpJpZesp7/MyZwt+XSiOIHt8MgOT1r6F7j1ENGB317/Vt3545FxaP+4cakKhEnG3fXDxuvcF8VoWvXLF5uljMFiRGR204eeuWtWnHjuyblzs6bPh5+nTq4w89VNz2lZWV7CXCcu0VpHAEEXzvAtglIuOwwPsigC8F8L5UAYaGbPvM5GRiAWFgwObs0tXV2YrzrVuZz/X22j9bWvw3Bovk3lPnmZuzG8s8a2vA739vizje3Obx4za8zHanxp49FtyXL9vn7tiRvcebr4aG7Kvw3A9YeYoOPlVdF5GvA/gDbDvLf6rqsU2+jarIrl22YLC4mHubCGA9p3feSR2C1tUBDz5o/37//VZCKn0+bceO/M7Ynj2b+vXsrM2leYsodXUWPocPA88/n/19OjoKL+DgHY1L3rgdidjvUXpPEih+wYSCF8g+PlX9LYDfBvFeVJnq6vK7lay311ZZp6YSix/335/o9fT0WB28Eydsq4y3fSbfhY30klTeOVtVG057Q9ylJet5Bllo9M4d4MMPLWxFbEvQ/v2JwN692x6fmkqU0Nq9mwVIKxFPblDgOjpyr7IODNivlRULqnzOBa+t2es6OlLDL3mztt9iS1C8azK9eTxV67kuLgKf+IQFnogF3a5d1issR1l7KgyDj0KTT72+uTmbs1tYsHBrb7fQ8c4C9/RYT6y/PzVAm5qCvTf30iX/xYtYzAIweUVYhKFX6Rh8FKjVVStFf/u2DW9HRws/ubG4mDpfuLFhw9fmZtu7d/Omzav19KSuMEejdmF5kIUScl1jmc8Vl1RZGHwUmFjMNicn94ympuyESD69L+9OjatXbbFgYcF/uLq8bIVHvVLzqnZszTsJMjRUePXkuTkr9rC6aoG6fbsNc3MdYytllWsqDQYf5eSdbfXO/w4N2cqt3z65kyczh4PeyYxnn839OarAu++mXkV58qQFmN8RvFgsEXwiiXnDrYjFLOSWlux4mWpiVTYet/YsLtrCSzRqbdu2LfWYXlcXT2VUIwYfZbWyYhP63ubj9XULijt3rMeVzu/+XMCGpGtruee9rl3L/P6WFuv99fenDpe9fYVHj1owjo9v/Vjf3FzqMHpmxgJ69257z5mZRGn9q1dt5bmnx17f0GDht2OHbdPhpUXVh8FHWV24kAi9ZHNzFmbpw1e/unaAhcRmFwHNzmY+1t9vjy8sJIJvackC0vv67l1ry8rK1godHDuWOoxeWrIh9MyM7btLDmHv5EVDg/36zGfy/xyqTKzATFktLm7tuaEh/9du3775lhW/3mBjow2rve/3Skb53YNx+nT+21eWlzOPkXlDd29eMfnonYj9JfDee9ZLPHrU/y8Eqh4MPsoq1z21fs/t3p1Z7bm7G9i3b/PPGh72HzJ2dgIvvgh89rPACy/Y136vW1vLf3U1Gs18D2+lOBq13mnywsnt2zbcXVuzYfD0NPD228HuE6Ty4lCXshoZsVJL6edPe3r8iwlEo3bfxsJCYjtLvqXoW1vtcqCjRxOf19SUWUG6udk/4CKR/O/xra9PnNVNtnNnoic3NGRVVVpbE71D71gcYI9duVK6kmFUWgw+yqqpCXj6adtAPDtrwTY4uPnxss5O+7VVg4O2Mjs3Z0F2332Zc4Pj4/73+Q4Obq1y8v79NuRNrjwzPm77/+7csTnDF1+06tMnT1rgelWrPQsLDL5qxeCjnDo6bB9ePJ44llVK0aj/5UgerzLMqVMWTpGI9cLyGU4na2gAnnnGNkR71ZK94W3yvrzHH8++MMOqK9WLwUd52WxVtpxGR21O8O7dRDmoQnV15S5ksH279fjSb7ZrbGRvr5ox+Khi3Lxpx92Wl21ucGws+/A1EilPjysatT2LH3yQGGLfd58NlbMVO6XKx/90VBFmZmy7iFdtZXbWtpA8+2zhx8+C0tpq4ectfPAWturH4KMtmZ211c7FRZsT27kz95xcPuJx21CcXpz07l0rPLp3b3HvHxQGXu2ooJkbqnTXrtkG3hs3rPczP2+3iGW7sjFfi4vZ76vwO9FBVCwGH+Xt1Cn/KxhPnSrufXMtTrCXRaXA4KO8+V0SBNi+t2y3q+WjpSVxGVG64WH/x4mKweCjnJaXba/bxkb2W9AaG/MrH5/LY4+lbiuJRKxYAIOPSoGLG+RrbQ04ciRxI1p9ffbV1a1URcmmqckuKVpYsLDt6sr/CBrRVjH4yNeRI3YW1bO2Zr96emwxwquvNz5uK7tBKfS4G9FWMPgow/Ky9fT8RCLApz9tq7CNjZV1oqNSra/bivjGhlVwDntfIjH4qt7GBnDxYuK+iZGR4u+AWFnxX70FLBQjkcIvEHLN9evAoUOJAq0iVuQhiOkBKhyDr4p5d70mF9WcnraFgmLOkba32zA2vRwVkH+ZKbLfv+TQA+wvlGPH7Ngbh/Th4UClik1PZ1YSVgU+/LC4IpmRCLBnT+bj9fXBzufVumvX/EvxA3abHIWHPb4qlu1yn9VV24JSTO9sfNyGs9PTNrzt7rbQYymm/OXa21jMvkcqHoOviuXaOxdE5ZBCrmxMF4vZ2d75eZuDHB3NfjdHrdm2zeb0/OZLiz3fTMVh8FWx4WH/s6wdHZVxyfXSklUw9qqaxGIWgEtLdj9HrWtutimDkydTH9++3UKRwsPgq2KDg1bDbno60atoabGqwZXg7Fn/28gmJ20oXUwB0Wqxa5ctZMzM2Lxrf7/94l284WLwVbl9+2xrxPy8Hejv66ucP1Q3b/o/vrFhm6BdWSHOdjkThYfBVwNaWrKfow1Tc3P2wgbcxEth4nYWKpmxMf/H028rIyo3Bh+VTF+fXdfoFRsQsYn9xx4Lt11EHOpSSQ0P2yJMLGZzkKy4QpWAwUclF4nYMbhyuXzZ9g7euZO4F4RXQVIyBh/VlJkZ4PDhxNe3b9vXqu5snKbNcY6PasqZM/6Pnz5d3nZQZWPwUc1Qtf2BfmKx4go3UG1h8FHNEMleRKG5mUVTKYH/K1BNeeCBrT1ObuLiBtWU0VEb8p45Y+W0mpos9MbHw24ZVRIGH9WcsTELwPV1K89VKWeXqXIw+KgmibhR/YUKwzk+InJOUcEnIt8TkZMi8oGI/EpEuoJqGBFRqRTb43sdwD5VPQDgNIDvFN8kIqLSKir4VPWPqurdI/VXADwUREQVL8g5vq8B+F22J0XkZRE5KCIHZ/0uiiAiKpNNV3VF5E8A/O7aekVVf33vNa8AWAfwWrb3UdVXAbwKABMTEz73ThERlcemwaeqL+R6XkS+CuAlAM+r+l2kR0RUWYraxyciLwL4NoDnVHUpmCYREZVWsXN8PwLQDuB1EXlfRH4SQJuIMty9C6ythd0KqhVF9fhUdWdQDSHyMzcHfPihlZsSsYuKDhxgCXsqDk9uUMWKxYB33knU2FMFrl4F/va3cNtF1Y/BRxVretq/eOitW3aBOlGhGHxUsZZyLJfdvVu+dlDtYfBRxerKcfK7s7N87aDaw+CjijU6aoVE0w0O2rWRRIViPT6qWA0NwDPP2A1ps7NWVHRoiGXkqXgMPqpoLS3Ao4+G3QqqNRzqEpFzGHxE5BwGHxE5h8FHRM5h8BGRcxh8ROQcBh8ROYfBR0TOYfARkXMYfETkHAYfVb3VVeDmTfsnUT54Vpeqlipw7Bhw/rwVLI1EgJERYN8+K1NPlA17fFS1zpwBpqYSVZrjcavafOpUqM2iKsDgo6o1Pe3/+PnzZW0GVSEGH1UlVWBlxf+51VX/uzqIPAw+qkoiQHe3/3NdXTbfR5QN//egqvXgg5kBF4nY40S5cFWXqlZvr5WmP3sWuHPH7uF44IHclxQRAQw+qnJdXcDjj4fdCqo2HOoSkXMYfETkHAYfETmHwUdEzmHwEZFzGHxE5BwGHxE5h8FHRM5h8BGRcxh8ROQcBh8ROYfBR0TOYfARkXMYfETkHAYfETmHwUdEzmHwEZFzGHxE5BwGHxE5J5DgE5FviYiKSG8Q70dEVEpFB5+IDAP4NIALxTeHiKj0gujx/QDAtwFoAO9FRFRyRQWfiHwBwIyqHgmoPUREJbfpvboi8icAAz5PvQLgu7Bh7qZE5GUALwPAyMjIFppIRBQsUS1shCoi+wH8L4Clew8NAbgM4AlVvZrreycmJvTgwYMFfS4RUTYickhVJzZ73aY9vmxU9SiAbUkfOA1gQlXnCn1PIqJy4D4+InJOwT2+dKo6FtR7ERGVEnt8ROQcBh8ROYfBR0TOYfARkXMYfETkHAYfETmHwUdEzmHwEZFzGHxE5BwGHxE5h8FHRM5h8BGRcxh8ROQcBh8ROYfBR0TOYfARkXMYfETkHAYfETmHwUdEzmHwEZFzGHxE5BwGHxE5R1S1/B8qMgvgfBk/shdALV90Xss/Xy3/bAB/vqCNqmrfZi8KJfjKTUQOqupE2O0olVr++Wr5ZwP484WFQ10icg6Dj4ic40rwvRp2A0qsln++Wv7ZAP58oXBijo+IKJkrPT4ior9j8BGRc5wLPhH5loioiPSG3ZagiMj3ROSkiHwgIr8Ska6w2xQEEXlRRE6JyKSI/GvY7QmSiAyLyJ9F5LiIHBORb4TdpqCJSFRE3hOR34TdlnROBZ+IDAP4NIALYbclYK8D2KeqBwCcBvCdkNtTNBGJAvgxgH8E8DCAfxKRh8NtVaDWAXxLVR8G8BSAf66xnw8AvgHgRNiN8ONU8AH4AYBvA6ipFR1V/aOqrt/78q8AhsJsT0CeADCpqudUdRXAzwF8IeQ2BUZVr6jq4Xv/vggLiMFwWxUcERkC8DkAPw27LX6cCT4R+QKAGVU9EnZbSuxrAH4XdiMCMAjgYtLXl1BDwZBMRMYAPAbgnXBbEqgfwjoZ8bAb4qcu7AYESUT+BGDA56lXAHwXNsytSrl+NlX99b3XvAIbQr1WzrZR4USkDcAvAHxTVW+H3Z4giMhLAK6r6iER+Yew2+OnpoJPVV/we1xE9gMYB3BERAAbCh4WkSdU9WoZm1iwbD+bR0S+CuAlAM9rbWzOnAEwnPT10L3HaoaI1MNC7zVV/WXY7QnQMwA+LyKfBdAEoENEfqaqXw65XX/n5AZmEZkGMKGqNVEVQ0ReBPB9AM+p6mzY7QmCiNTBFmqehwXeuwC+pKrHQm1YQMT+Bv4vAPOq+s2w21Mq93p8/6KqL4XdlmTOzPHVuB8BaAfwuoi8LyI/CbtBxbq3WPN1AH+ATfz/d62E3j3PAPgKgE/e+2/2/r0eEpWBkz0+InIbe3xE5BwGHxE5h8FHRM5h8BGRcxh8ROQcBh8ROYfBR0TO+X+sucJdUmbi0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calc and plot whitened coordinates\n",
    "B = np.diag(1/np.sqrt(L)).dot(A)\n",
    "plt.figure(figsize=(5,5)); plt.xlim(-5,5); plt.ylim(-5,5);\n",
    "plt.scatter(B[0,:],B[1,:], marker='o',color='b', s=50, alpha=0.3, edgecolor='none');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoYAAAE/CAYAAADbpwJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xlw3Ol93/nP0wAIAsRNAAQBggDBc26NhnPSI49Gki3bir2brXISbw472Z04u8natUqp1lbtliu1qfJWqhJvyt44inO4NtrYieOsj5LKGq0la3TMQY40vG+C4I2LIEDcQD/7x7d/6uvXjQa60ef7VYUS+9fdv36aI37xfa7v47z3AgAAACKlbgAAAADKA4khAAAAJJEYAgAAIIbEEAAAAJJIDAEAABBDYggAAABJJIYoMOecd84dyvDcf+uc+1qx2xTGOffvnHP/e6nbAQAbcc79tnPufy30aze4z3Asntfney9UFv6Do2i891+W9OXgsXPOSzrsvb9aulYBQHnz3v/idrwWCMOIIQAAZco5V1fqNqC2kBgiJ865X3DO/UnC4yvOuf+U8PiWc+5jsYefjj0/45z7Leeci73m551z3479+Vux137knHvsnPsrseufc879IPbe7zrnnk34jFHn3D90zp12zj1yzv2+c25nwvPZ3vu8c+5D59ycc+73Jf3wfQBQbM65J5xz34zFq3POuZ+OXf93zrl/4Zz7inNuXtInU5e+OOe+4Jy755y765z77xKX8CS+1jn3hnPutnPu88658dh7fiHhPj/lnPu+c242FsN/rbh/CyhHJIbI1V9Iet05F3HO9UvaIelVSXLOjUhqkXQ69trPSXpR0rOSflbSj6fezHv/idgfn/Pet3jvf98597ykfyPp70raLelfSvpj51xjwlt/VtJnJR2I3f/nY23I+F7n3A5J/6+k/1tSl6T/JOm/yfcvBAC2wjnXIOlPJH1NUq+kfyDpy865o7GX/JykfyypVdK3U977WUn/s6RPSzok6Y0NPq5PUrukAUl/R9JvOec6Y8/NS/qbkjok/ZSkv+ec+6/y+W6ofCSGyIn3/rqkOUkfk/QJSX8m6a5z7pikH5X0jvc+Gnv5r3vvZ7z3Y5K+EXtPLt6S9C+99+9579e9978raVnSKwmv+efe+7ve+2lZYP1YDu99RVKDpN/w3q967/9A0gdb+osAgPy9IutM/7r3fsV7/+eS/lTSX4s9/0fe++9476Pe+6WU9/6spH/rvT/nvV+Q9GsbfNaqpH8Ui31fkfRY0lFJ8t5/03t/JvY5pyX9B1k8Rw1j8wk24y9kvdNDsT/PyILIq7HHgfsJf16QBcBcDEn6W865f5BwbYek/iz3Dp7L9l4v6Y733ic8dzPHNgFAofVLupXQmZYsJg3E/nxrg/eeTHic7bWSNOW9X0t4/MOY7Jx7WdKvS3paFi8bZTMqqGGMGGIzgsTw9dif/0KWGP6okhPDrbol6R977zsSfpq99/8hz/fekzQQrHWM2V+A9gLAVtyVNOicS/wdvF/SndifffpbfuiepH0JjwfzaMf/I+mPJQ1679sl/bYkl/0tqHYkhtiMv5D0SUlN3vvbkt6RrffbLen7W7jfA0kjCY//laRfdM697Myu2OLo1hzule2935O0Jul/cs41OOf+sqSXttBeACiE92Qjd1+IxaQ3JP0lSb+Xw3v/o6RfiG1eaZaUT83CVknT3vsl59xLsrWNqHEkhsiZ9/6ybH3KO7HHs5KuS/qO9359C7f8NUm/G9uV97Pe+5OS/ntJvynpoaSrim0uyaFtGd/rvV+R9Jdjj6cl/RVJf7iF9gJA3mIx6S9J+glJk5L+L0l/03t/MYf3flXSP5et374q6d3YU8tbaMr/IOkfOefmJP1vsqQTNc4lL7sCAACVwjn3hKSzkhpT1hICW8KIIQAAFcQ591/HSnF1Svo/JP0JSSEKhcQQAIDK8ncljUu6Jmld0t8rbXNQTZhKBgAAgCRGDAEAABBDYggAAABJJTr5pLu72w8PD5fiowFUuFOnTk1673tK3Y5SIHYC2KpcY2dJEsPh4WGdPHly4xcCQArnXM0eZ0jsBLBVucZOppIBAAAgicQQAAAAMSSGAAAAkERiCAAAgBgSQwAAAEgiMQQAAEAMiSEAAAAkkRgCAAAghsQQAAAAkkgMAQAAEENiCAAAAEkkhgAAAIghMQQAAIAkEkMAAADEkBgCAABAEokhAAAAYkgMAQAAIInEEAAAADEkhgAAAJBEYggAAIAYEkMAAABIKmBi6Jyrc8593zn3p4W6JwBUO2IngHJSyBHDX5J0oYD3A4BaQOwEUDYKkhg65/ZJ+ilJv1OI+wFALSB2Aig3hRox/A1JX5AULdD9AKAWEDsBlJW8E0Pn3OckjXvvT23wureccyedcycnJiby/VgAqGjETgDlqBAjhick/bRzblTS70l60zn371Nf5L3/kvf+uPf+eE9PTwE+FgAqGrETQNnJOzH03v+K936f935Y0l+V9Ofe+7+ed8sAoIoROwGUI+oYAgAAQJJUX8ibee+/KembhbwnAFQ7YieAcsGIIQAAACSRGAIAACCGxBAAAACSSAwBAAAQQ2IIAAAASSSGAAAAiCExBAAAgCQSQwAAAMSQGAIAAEASiSEAAABiSAwBAAAgicQQAAAAMSSGAAAAkERiCAAAgBgSQwAAAEgiMQQAAEAMiSEAAAAkkRgCAAAghsQQAAAAkkgMAQAAEENiCAAAAEkkhgAAAIghMQQAAIAkEkMAAADEkBgCAABAEokhAAAAYkgMAQAAIInEEAAAADEkhgAAAJBEYggAAIAYEkMAAABIIjEEAABADIkhAAAAJJEYAgAAIIbEEAAAAJJIDAEAABBDYggAAABJJIYAAACIITEEAACAJBJDAAAAxJAYAgAAQBKJIQAAAGJIDAEAACCJxBAAAAAxJIYAAACQRGIIAACAmLwTQ+fcoHPuG8658865c865XypEwwCgmhE7AZSj+gLcY03S5733HzrnWiWdcs697b0/X4B7AyURjUpra1JDg+RcqVuDKkXsBFB28k4Mvff3JN2L/XnOOXdB0oAkghsqjvfSlSvS9evS6qq0c6d05Ig0NFTqlqHaEDsBlKNCjBj+kHNuWNLzkt4r5H2BYrlyRbp0Kf54aUk6fVqqr5cGBkrXLlQ3YieAclGwzSfOuRZJ/1nSL3vvZ0Oef8s5d9I5d3JiYqJQHwsUTDRqI4Vhrl0rbltQO4idAMpJQRJD51yDLLB92Xv/h2Gv8d5/yXt/3Ht/vKenpxAfCxTU6qr9hJmfL25bUBuInQDKTSF2JTtJ/1rSBe/9P82/SahV6+vS3Fzm5Gy77dhhawrDtLUVty2ofsROAOWoEGsMT0j6G5LOOOd+ELv2q977rxTg3qgRly/bdO3amhSJ2GaPJ5+0PxeLc9Lhw9KZM+nXjxwpXjtQM4idAMpOIXYlf1sSBT2wZaOjyRs+olHpxg2prk564onitmV42DaaXLsmLSzYSOGRIxIzeCg0YieAclTQXcnARlZXpVu3pMVFqb1d6u+3JDDMzZvS0aPFHTWUpH377AcAgFpDYoiiefRIevddaWUlfu3qVWl2Njz5W121dYfFTgwBAKhVJIYomjNnkpNCyTabzM3Z6GGqXbvs5JFEi4vS2Fh8xHFw0KZ+AQDlZ21NWl6Wmpro5FcKfqWiKJaWpIcPw58LAkY0mnz96NHkx9PTNuK4vm6Pb92yaegTJ6TGxsK3GQAq2ePH0uSkdbD7+mzddrFEo9LZsxano1Gr+nDkiHTgQPHagK0hMURRZDtvuKND+vjH49PKu3ZJBw+mb/g4cyaeFAbm521H8zPPFL7NAFCpzp5NXr+9Y4f00ktSZ2fxPv/mzfjjlRW71thoa8tRvkgMURSNjVJ3t/VeU/X3W7B68cXM719ctKQxzIMHJIYAalM0astr7t61s977+y0JTN3Ut7IinTolfepT2TvqhbC2ZiOFYW7cIDEsdySGKJpnn5W+9z1L8gJdXVY7cCPZ1qYUc3oEAMrJyZPWOQ5MT0sTE+ElthYXpakp66Rvp+Xl9KVBgYWF7f1s5I/EEEWza5f05pvSvXsWHDo6cq8P2Nhorw07KnZgoLDtBIBKMDmZnBQGpqdtXWFHR/pzmRK2Qmpqspi9vJz+XFibUF7YI4SiikQskTt8ePNFo597TmppSb7W1ycdOlS49gFApZieDr/e0WHVHlI1NEi7dydfe/zY1ndfv26bBAshEgk/LSoSyW2GCKXFiCEqRlOT9MYbNmoYlKuh9wmgVu3YEX69pye9NJhz0tNPJy+9uXw5+dSp8+etAz44mH/bhofjax2DGaLDh4nZlYDEEBXFOam3t9StAIDSGxiQLlywzR6JGhqkH//x+HrDHTvsNKfW1vhrZmaSk0LJNq+cPm0xthAlwPr72WhSiZhKBgCgAjU0SC+/LDU3x681NVlZmp07LSl77jk7cz4xKZRsF3OYaNTWgaN2MWIIAECF6uqyTX2PHtmIX0dHbuVovN/ac6h+jBiibKyuWk91fLw4O+cAoBo4ZwlhZ2fuNQr37s18r76+wrUNlYcRQ5SFsTGrih+cbNLYKL3wQvoOOgBA/rq67Hi61ELYTz5p09GoXSSGKLnZWVvwnDh9sbwsffCB9JnPlGcB66Ulq+y/tGQBdu9eDogHUFmefto2sNy/byOFAwPpaxFRe0gMUTATE/EA099vCVMubt+2pNB7K2vgnC2mXl21+5VbAevJSen99+Ojm6OjNo3z6qtSPf+iAFSQzs7inZ+MysCvMRTE6dPJB6ZfuWIjfQ0NlkD19toURWqBaskSwEePbEojqL3V3CyNjNhz5cR76aOP4klhYGbGCsSGFXUFgFKYmIjvPt67l1JfyA2JITbl+nUbIQumUI8etY0iiUmhZIVTp6akp56y1z14YMnTG2+kF2VtabHXJ04lLyxYja2f/Mnt/kabMzeX+azP+/dJDAFs3eKiLUkJaghGo9K1azarsr4u7dljMSaXGoNnzyavHxwbk4aG7Mx6IBsSQ+Ts4kUbCQxMTIQfyP7ggQWllRVpft56qgcOSG1tFpxSj7BbW7N1LbOzydd377ZErJzWvGRbR5jrbkAASDQ9LZ05E4+BPT1Wf/DsWetwBq5ckd57z5bq1NXZ/x49mt7ZDmZgUt28aaeaMHWMbEgMkZPVVRstTBWNWm82CEzr69K5c/Ep4UjENpJcuWI91bDzO5eXrRc8NWWjipGIJYUdHYU7u7NQWloswU1NYqXyWwsJoPwtLVmyl3h6ycSE9PWv2yxK0OGMRm3JzsSElfQ6csRmb6anpU98Irlj+uCB3e/uXYurks3c9PfbcySGyIbEEDlZWEhfVxfYuTNed3B6Oj6q5pwlUvPzdlC7JO3fn/7+ri7ryfb02E+icixX8/zzFsgTk9a+PjsbFAA2Y2ws/Ug7yZK/HTviZwtfvGjH30nW8V5aisfT+/eT6xJGIvb6xGUvDx5Yx/ypp7bne6B6kBgiJ01NFmzCCk8HCd25czay2Nxsr+/osN5tMErY0GBBsL3dDlMP9PfbtMfMTPJ9+/vtteWmrU361KesGPfyMrv6AGzd4mL49R07LL5IluCNjsafq6+30cSxMYuzs7PJiWF9ffh9FxepnICN8X8R5CQ4hH1sLPm6c9LBgzayt3ev9VIbGqRnnrHp42B9zK5ddq2hwTaV7NsXL6IaiVipl+vX7fV1dTYtOzRU3O+4GZEIU8cA8pep85t4fXraYqdkMXfHjngd1ceP02PR8rJVdbhxI96Zd85mNYJlPuVgcdES3rk5+x0xPGz/i9IiMUTOnnnGepvB1EdLi3TsWHy6d+dO6WMfs8Dz4IFtGgk2mrS22gigZD3d8fHkxK++3tbMsKsXQC3Zt886xfPzydf7+izmfvSRreNubbUksKnJ4mtQxSEatcdXrsRnYpqbLS63t8dnYjo6LM42Nxfvu2UzNyd95zvJJclu3rRBAmZgSovEEDmLRGx9yhNPWGKYuhMucPy49QIfPrRg19VlZRYSF0dzSggAWLJ24oSV7ApmTPbts051JCK98oqN/n3rW9bxfuedeFLY3GxHh9bVWVmbkZH4buWLF21EMbFqRGNj+cx0XLiQXqd2fV06f97+PlA6JIbYtEgkc1IYPD8yYgnhO++kP19XxyHtABBobLTRwWeeCX++t9eOr7tyxdZzBwcFDA/HN6esrtrUbEuLxdjXXrNdzJOT9vzu3VYZolzWGE5MhF+fnrYEsRyPQq0VZfJ/EVSjjg477eTChXgPt75e+vjH4+tlAAAbO3bMRgLn5216uLMzeVq4rs6W8wR27bJp2WBNYbbOfCnU14evd6yroyZsqZEYYlsdPGhTF+PjNpLY11c+PVYA2Mjjx5bAtLeXfhSrrU16800rl5VqaCg8tpZbQhjYv1+6ejX9+sAAS41KjV/R2HY7d4bXL6xmc3PxXyYkwkDlWVqSTp2yqU3J/h0fO2anOJVSb6/Nuly6ZKOHDQ2WFB47Vtp2bdbRo9b+e/fi13p6qLNYDviVBRTQ4qL9Mnn40B6Xyy8TAJuTmBRKtuHu7Flbw5daiL/YBgbsZ2XFYkwljrBFIrZR8fFjq8MYnCqF0iMxxLYI6lPNz1uZhaGh5PUv1erkyeRC3cEvk9bW9DOlAZSnubnkpDBRcEpTKSwtWemapSXb3Ld3b+Wvx2tpiW+mQXkgMaxw0agdmF5fb8lHOZiZkb73vfgxT/fuWZJ44kTxAsDKigXwhw+t7tfw8Pb//Tx6lH56S+DmTRJDoFIEJ45s9rntNDkpvf9+/GjSGzcsOXzlldKvfUR1ITGsYPfuSWfOxANVW5vVtCp17+vcufSzP1dWbHfyiy9u/+cvLkrf/nbyWcZjYzZtsWfP9n1uthMFSvXLBMDmBRtNws6HL8X57d5boevU9kxPW3HsxCNGgXxV4MoESDbVcepUcsIxO2u71YLSMKWwvp55CmZ8fOv3jUbtvplG5BJduZKcFAbvP3du65+fi46OzD33UvwyAbA1DQ3hpzA1NZVmvfDsrJ2XHCY4dhQoFEYMK9TYWHgCuLBghUN7e4vfJsnWu0Qi8fM5E211d+69e9ZbfvzY7tHebiOjmRYqZyqcOj9vP9t1FmdDg/XcL15Mvt7czOYToBxNT1vnurMzfQ30oUM2+zI2ZktSIhH7d1yKadts6wiLscYwGrUEdGnJ/q44sq66kRhWqHJcAyNZ8BwYsAPeUw0Obv5+jx9LX/2qLbheXbX79/TYtO1nPhO+Gy9bArrdpWMOH7a1jGNj1sbubjsFplxriQG1aHHR1uvNztpj5+zf6ZNPJr+ur8/W9j14YI/Pn7cZiRdf3J5ZgLU1S8DW1qxzHxSwbmuzuDI3l/6e4Az6zVhashi1sGAd7cHBzLFxbk56993kWZi+Puucb+du6Olp+3sPfqeUeolULSExrFBdXdKdO5mfK6WnnrKAMzUVv7Znj9Wt2qxTp2yRdSA4MN576bnnwo/WGxwMnzbu7bWjp7ZbXx9H/gHl7MMP40mhZPHk2jVLwPbti1+/fz85/kjWQT11Svr0pwubGE1MWFWDxPXZhw/H6xM+/7wlaIlrmfv6bGPdZjx8KL39tn239fV4Yvj66+GVI37wg/SlOcHfy8GDm/vsXJ05YxsWA5cv25GAzLwUB4lhhRoctJ2uicFNsrIw2zVVGlhdtc+enLTp0/37k8s3NDTYOZ0zMzbi19a29fpUqUE5MDFh9w5z4ID1cm/dik+3d3TYAfQAatv8fOZ10GNjyYlhps738rLFv9QlO9GozW4kjnTl0klcX7dkM3XT3pUrNuvQ3W0J3Kc/bUtrgnI1WxkE+NrXrIRWYHravkt3tyWfiYLj98LcubM9ieHkZHJSGDh3zsrz1ELZs1IjMaxQwSHpN25YEKqvt4C2lenazVhdlb7zneQpjbt3bQomNUh0dMQPeN+qTFOw3meeWnDORhMPH7ag1tTEmhgAZnU19+fC1kpnes57m55OXON89651VJ9+Onubxsczt+vOnXipq7q65MR1s+bmwmdTguupieFmvn+hZNpM4739rhsa2p7PRRyJYQULds6NjFjwePjQerKDg9s3ZXrzZvg6l0uXbOSwoaGwn3fwoPXi5+eTr/f3b7zBprk5+ZD5wOystXdqyhLP/fvtcyq9UCyAjbW1WXwMW4udWri6ry88UamvT69Lev9++Ma3GzdsujfbGrlsSVZYyZytyjRSKlkd1lStrTYDlRp/pdIslyFGFwflasrY6qr9g8wWNJaWpG99Szp92pK2Cxekb3wj/B95IUxOhl/PVqYmH0eP2mjkvn3xI5NGRmzjSbYSDpksLEjf/a4F8eDv98IFW9MCoPpFIrYOOjXJaG5On/UYGEjvgDonPfNM+maNTNUQNnpOsoQ003rFQiZgTU02JR0m0yaW555L34nd3r596wsztSMSYe12sTBiWIbW1y1RuXPHksLGRkuQwobQg4PUE62u2vt/5EfC7z89LV29amv0WlvtH3iua1WyjQhux87bjg5bFH3lSnxaeMcOWyh+4YK9pqfHpkByGSW9fj18ymZszEZfWb8CVL+BARsJu3kzvl5veDg9vkUi0ksv2RTmxIQ9H3RSU2WLjanPLS/bzEswq7FjhyWrqR3Uvj5bV1co3d22meWjj2xndqC93Q4AWFlJj+O7d0tvvmlrtoNyNf3927cjuavLlgFduRK/FolYgkp1h+IgMSxDp0/bAubA8rJda2xM7zEFZRRSBdPKqcnS+Litgwk2ZczP2z1efjm38z/377d1M6laW7dvHV9npwVnyUYsv/e95OcnJmzh9muvbXyvTCOp3lugJjEEakOua6Cdy63SwOCgdVhT68s2NMTf672t5bt5Mz4TtHevdWyHhy0pun07Xq5mz57CTp9GIhYn6+vtd8Hyso0GNjXZ7mPJksdnn03exLhzZ/jpKouL1kG/f9/u3d8vPfFE/kuKjh2zBPz+fWsfm06Ki6nkMrOyEp54SeE7dDP12oJC06kuXkwPXN6nF2XOpKfHpnYTpxZaW4tz1J1kATXM1FT42sdU2XZsh61HBIBctLRY5YPEKebGRouNwbXr1y2OJy4Puncvvku4rc3i67PPWjK5HWvqOjttKc6P/Zjtcu7pSY6LQed7o7WNa2u2EfHOHXttUK3i3XcL086WFisyfuAASWGxFWTE0Dn3WUn/p6Q6Sb/jvf/1Qty3Fi0tZV5TmDj0HxgYsGnhVL296b22aDTziNnMjCWIuQSigwdt5HB62ob2i7njN7WeVqLlZUtSwwQlGYJTWVKT5r6+7S/zA6QidlaXffsslkxNWYzZvTs51mTq2N6+bTuXU9fyTU0lT3cfOFCYjYWRiI1Gjo6Gn6C1uGgJa7Yd0Ldvh/9OmpmxWZxcZqBQnvJODJ1zdZJ+S9JnJN2W9IFz7o+99+fzvXct2rXLErqwdXBh0x5Hjtg/xMRNIS0ttjg6VSSSeTfezp3JSeHysk2LJNYqHBiIP9/QYIGl2Lq6wje51NWFL6r23orZJo7CLizY6GAkEq81tlE5CaDQiJ3Vqb4+c2xMLE6dKBq1EbjExHBszNYCBqambJ1fUIg6GrW4Fqx9HBzMvLEkk2yb9zba2JdthmZ2lsSwkhVixPAlSVe999clyTn3e5J+RhLBbQvq6mwtx/nz6dcPHQp//auvWrI0M2OJZW9v5pG/4WHbsBJ2PbCyIn3728mBYXLS/rE/8UTmts/OWs+2o2P7FgkfOGA91dSRw0OHwte13L6dPjXf0mKB9fXX7T2lOPsUELGz5uzeHV7+pqUleSQwGo1vrku0tGQd9ieekN57L3lA4MYNm4LeTJ2/bIlkpucePrQEdXTUPj91VFTi+LpKV4jEcEBS4sm4tyW9XID71qyDBy1xGR21ofrOTksWs50ekmsV/MOH42tB1tctKRoeTk46R0fDe4vXr1upmNSpjKUlO8rp4UN7HInY/RKPwBsdTf4+R49ubQo6SOiuXo3XIRwaylziINN6zaUl23iTz3mnk5OWZAe7pUdGNn88FWoasbPGHD1qcSPxhBPn0s9onp3NPLo4MWFLZsJKh507Z7Ew180fe/fa75XUE7Q6OsLrxN64EV8Pub5uHe/JSfteQXLY2rpxjVmUt6LtSnbOvSXpLUnav39/sT62Yg0MJE/dFopzVhbhyBFL0pqb0+txZapHGI1a8pe6O+/DD+NJYfC6y5ctQPT3258TRyknJiyp+5Ef2fzUh2TJYa5Tv2HrZ3J5biPT07bIOnF395kzlnSH7d4DtorYWT3a2qRPfMI62Y8eWfwdGUlfJpRtxqWxMXM1ivV1i6251vuLRGzG6fJlW1PonCWLR46kzzqtribPZNXX2+7hmzftM/fssc8NqxGJylKIxPCOpMSD2PbFriXx3n9J0pck6fjx43n8SkYhNDRk7lVm2wGW+tz8vAWFMGNj1nO8di39uWjURv1eeCG39koWmCKRzU399veHF5dtbNzaOaOBK1fCE8tr12zEd7tqfKGqEDtr0K5d4WvAEzU32xq9sNi1f3/mxFDKLT4uLcWX5HR12YjlRp3tycn0jZFNTZYc9vRYSTHiXnUoRGL4gaTDzrkDsqD2VyX9XAHuixIZGrI1JKna29N7tpmmOyTbwLKwkH4wfCDX01kePbLpi+npePX7oSGbwq2vt+QvUw87qIWVGEiDYqn5BLHUqZfA6qqNxLLDGTkgdiKj55+3JTrBDE6wRGdgwDr1d9K6ENZxTz2qL9XkpNWyDcrRBGXQ+vstnu7dGz4VnDqzlKixMb94urJiBy40N1OaphzknRh679ecc39f0p/JSi78G+99yDHdqBSdnVaP6/z5eOLX1SV9/OPpr21vt6QsLEHs6bEeZVAiJlUuydPSktXUCnZpR6N2pN3bb9uUhWTtfOGF8J2AwckFExPxHdb79uUffHbtCi+dU1e3fedUo7oQO5FNY6N04oR1QpeX47FWssTt8GGbdQlmLhob7fSSjaZxP/oouUbh6KgVu963z5LDsTHreD/7bPL7urstbobFvWxlbbLx3uL36KhRyZQrAAAYgUlEQVTFduesDWHH8KF4CrLG0Hv/FUlfKcS9UB4GB61n+uiRBaNMSVwkYjvkEssqSPHNGEEZhbD6XSMjG7djbCy5dM/DhxbEJAuYbW0W5L7/fSvamimY9PRsrXyC9zbimTrtfvBg+BT68HD2njWQiNhZe6JRi10NDdk7x9Gorfubn7c4lxqDjh2zBC4oV7Nnz8ajdqnnyz9+HI+n09PxTXw3b9qUdeIMkXNWrPv99+Mlz5yz9YhbLU1z/br9BLy3kdCGho2n27F9+BWGjCKR3HYO799vAW501AJGaiHWoHDr2JglWbt22S62XIJJ6jnQiRtjEusxrq5agCzUIeve24LsGzfs3kGbgw1Be/bYCGpwVnVDgyWFiTuxASDR7ds2QhbErt27LY6kzmAsLNhMSWIS195uG0USE8SmJou/2ayvW+ydmLDPffQovulvZib+utSk8sGD9KVDHR12Wsr4uMXyYBRxqzIV/L51y2aEWLNYGiSGKIjduzOXfolE7B/5E09YkpXLVOv6uvWYs5XoaWpKfpzpxJituHgx+USZ+XnbfZ1YvDbYOb66atfZiQcgk4cP7TzixE1rU1PSBx9YCa5EZ86klwx79MjiUi4jaWtrlgTu2GH1DhOrRty6ZfGsvz95hiV1M16m2ZdgnXchhB22IFn8X1vbvnq4yI7EEEUTnLySzfKybTS5f98SvdbWeL1FyXqsU1OWMCYWUa2rK1yl/fV1G/0Mc+1a+lrGfA+MB1D9Mh0/NzNjP8Ho3Npa+G5kyeqyZksMo1GrZTg2Zn+enraYundv/DUjIzYb0tNjyeDt2+l1C50rXLm0+Xkr1j0xYR3owUGbfg6ODAzbYd3SQlJYSiSGyNvCggWiYC3M0NDW/1G/917ybuXg2KXubguee/fGjw0MOGcLpQuVoC0vZ95JnTq1DQC5yDQ6JqVv6MhUY3Wj2qsXLiR3aqembLSwvj7ecW5utuRyYMASw2PHLH4HG1Lq6mzzR+qMzFYsL9spWsHmxLU1K/U1N2frFYOC34mbYZzLfsIWth+JIfIyPW3JXJBI3b1rgenECQtAG5mdjZ+04r31KlOTyro661m++mr82sSErXNpaLAAt5XyMOvrtsblwYP4mcm7dllPu74+PDlsbd385wBAV1f4SGDqWu4giQt7beLIX6ogniUKNsI9eJA8o7K6avFtft463UeOxD+vt7dwnezR0fCKFffvW3LY3m4Fv69ds98FQcHvrZyKhcIhMUSa5WVbIB0cJ9fXZwVQw3qQ586lJ1BLS7YWJqy8TaLxcVtfE6wNHB+3HWlPPBFeSDvRVnYZB4fQr65aMLx9O774em7OSuC0tdnGmakpC6qJJ7M4x6kmALZmeNjiT+rawbBjRp9+2jafJI4ktrTY6F4mq6vJI2+Sxblg00kgOOu4sdFG627etI73yy/nXiJmbc1i9cKCxchdu2x6enw8Xlt2797shbhnZ62j3dJiI5QoHySGSOK9BaRgCleyBHFmRnrjjeTAsbKSvKstUVACIZtz55I3jDQ3W3C7c8fKwSR+zr17NjpZV2c1s0ZGNrdj7erV5EPpz52zoHXsWHwH8vp6fGfz7t3xs5jX1mwNzpEj+Z2tDKB27dhhx4BeuxafGRkcDK8B2NIiffKTFguDJTr9/dljXmNjep3B1lbbtfz4sT2ORi2eHzyYvFluasoSxFxKiM3N2e+IINlcWbHqDQcPWlJ4+7Y939ZmbVpYsOdSRyE5BKB8kRgiSTDEn2phwQLKYMIBXpGIBZewdS8b1fJbXIwHq0BLS/qB7mtrFkgPH473hi9csHWIuR6nt7ycfFazZJ/x+HG8FmFiT3t21kZHd++2gHboUG6fAwDZNDba7Esu6uttvXaunLM1e6k1Zfv7bT1fUD/RufAKCg8e5JYYnj6dPAI5Pm7x+O5d60AHaxwfPrTO9Pi4dcwT1w12daWXwkH5IDFEkrCkMNNzwfFJwZRzoo0q4QflXVKTysOHLaA0Nlqy5r0Fq9RE8+5de222cjaBsDM+g/vNzKT3XBcWbCp9cdECXktL4cozAMB22b/fRiKvX49P8x46FF+z19iY3kkO5DIDs7ycXEtWsjgpxXdAJ1pft8+/dSt+VGh//8bnMqO0SAyRJNvmisTyMIFnnrEAlDil3NdnPcVsGhrsdffuJV+vq7PFyMFU8ocfhp8JKtln5pIYhi2k7u620dG6uuTvvLgY34wiWUL5wQc2OhmcCgAAhTI7ax3gtrbC1ELt68vcke3osCU7qescpa2Xp0lcD554StXqqk0rr6zY92tqsqn0xHXbKE8khkiyZ48lgKnTvE1N4YFjxw4rzjo9bWth2ttzS9YkKzGzshKfznXOpqoTpzOylUzItZxC2BmfwZF9LS3Wi+7stGCZOIrZ1haf7rh8mcQQQOE8emQd3yDWNjfbJozu7vzum+1ceOdsU+D77yfvFs601jFVY2N8/XWgt9emi7u67P6PH8drKAad8kjElgW995705pscG1ru+M+DJJGIlYU5f95G87y3ZPGpp7LvWOvqSq+cv5EdO6TXXouf39nenp7sDQ3ZwubU3XatrbkH0EjE1th88EE8OYxELKE9cMAC2IkTFqj/4A/sO3d2WsALevBzcxbsOKIJQL7W16V3301OzhYWLGF7882tHTPnvXTqVPIszKVLNtuROILY2WnH2t2/b1O/3d25d+Yl69An7phubJReesmmiScnLZ5Go8nxee9eSwaXl20UcXh4898PxUNiiDQ7d1qvMhg52+6j3traMgem5mYLOmfOxHvWPT3Ws95Mu4IzPicn4+VqgnqJQQDr6pKefz6+ZibRzp0khQAK4+7d8Pp+6+uWOG1lw9vt2+lLc6JR24zS25scv+rqtj513NJiyevdu/EO/Z49Fo/X1qSf+Anpm9+0zw1qMiZWc+CQgPJHYoiMyuXs3+5uK92wsGABLZezlsM4t3Htw5ERK2UTdh0AturxYzthZHnZZkkSj/pMlO2ElGxSk8JAsFynUEeGStbuxAoVgfp6+3nhhfRZngBrDMsfiSFKZm3Namcl1vTKFrxyOUklXyMj1su+ds0C6o4ddi2xriKA6ra0ZDFgaso6okND+VUmuH/fpnmD6ghzcxb7jh1LX2+32SU5gXLpyEs2QtnRkV7ntqWFtdqVgMQQJbG2Jn3nO8k1C4NTT0pdN/DQIUsGg8SQKWSgdiwtSe+8k7xZbXzc6g9upYPovS2FSSyZ1dpqCef9+8mbPrq6tp6A9vfb/VIFG0aKyTnplVdsjePdu/Z30NdniTDxtPyRGGJD09O2K/fRI1tgPDKSf6/v5s3kpDBw6ZL1zgt1VudWRSJbWwAOoLJdu5acFAYuX7bYtNkdtY8ehd/v4EG73tkZT5xGRrY+8tffb6W2Est71dXZuulSJGMNDVavkJqFlYfEEFlNTdkOtGAjysqKTYmsrm6uKn+qycnw69GofSYFpQGUQmoB58DamnVmNzvVm6maQyRiy2dee21z98skKEVz4IAtz2losA0mwSY7IFckhsjq8uXwI+8uX7Yq+1vt3WYbESSQASiVbJvbtrLxrbXVNlw8epT+XC61AzerszN+0gmwFcz2I6uwYCbZFMhWd89JllSGaWnZ+uJrAMhXphp73d3px2fm6uMfT988t39/5jgIlBIjhkizumq1ppqb7ScsOayvz29kr7vbimZfvBgva9DSYoWoAaBUenttXdylS/Ej3np6bK3eVgW1/8bHrUPd1RV+xChQDkgM8UPe24kno6PxUz527Ag/8WN4OP8FzSMjtsYmODqJkUIA5eDAARvNm521GLjVkcJEzlkh6K16/Ng2lkSjtgab6WJsFxJD/NDVq9L16/HH0ahNGe/aZX9eXraRwqEhKztQCA0N+QVLANgOdXXlk3zduCGdPRt/fPWqJa/s+MV2IDHED42Ohl9fW5M+85l4Xb9sZybXuvV1my4Kjt0rRlFuANVraSn8NKYbN+wM4mLXKCwG720goqGB3zelQGJYoZaXrRbg3JytVRkayr/uXqbNJME6m6am/O5f7aanpQ8+iJ+B6pwVyy7U6CqA2nP/fnhliOC5aksMb9+2teeLi5YU7t9vxcUpjF08JIYVaG5O+u53kw9hv3FDevXV/M6h7OwMr+HV1kavbSPRqNV3TPxv4r105YqtneztLV3bAFSucjrqbrs9eCB9//vxx+vr9rstGpWefbZ07ao15OAV6MKF5AREslG98+fzu2/YcUXOMeKVi6mp8NMNJOsBA8BW9PVlHi2rtnOHb9wIv37rVnzmCtuPxLACjY+HX5+cTD6Pc7N275ZOnLBg09pqAem119gckoug5M5mnwOAbBobbbQsdeTwyJHy2RxTKPPz4deDjZAoDqaSK1B9fXjvqa4u/2mHjg7phRfyu8fkpJ03+vixJZiHDlV/KZrdu+3vPywJJLEGkI/BQaulePduvFxNNdZBbG+XFhbSrzc0sJGvmBgxrECZjlEaGCj9epT796V337VRzYUFWzPy3e/a2Z3VrKHBCnan6u7enmOvANSWnTut9uuhQ9WZFEr23cKmzQ8eZJ17MTFiWIGOHbPRuMRka/du27lVahcvpu+g895OEejpKU2bimVoyKZ2gvUwPT1WToLddACwsY4OW750+bL08KFVwgiKjaN4SAwrUH299Mor0sxMvFxNOaw1WV+39oSZmSluW0qlrS185BAAKsGtW9LYWLwW66FD+ZdC24zOTunll4v3eUhHYljBOjrsp1wER+il7piWihtYAACbd/68rQ8PzM3Z8qDXX7dNMKgNTHKhYJyzM5TDHDhQ1KYAADZheTm8XMziYuZTsVCdGDFEQR05Ykfo3bxpU8v19ZYUHjxY6pYBADKZmclc7uzhw+K2BaVFYoiCcs7W2B09aj3NpiZLDgEA5Svbch+WAtUWfmVjW9TXWw3DYrh3z46eCzbiHDpkpXsAALlpb7eNH6mjg85ZxQXUDtYYoqLdvSudPCk9emTTILOz0ocf2s46AEDujh9PLivW2Cg9/3x5VL1A8TBiiIp25Urm64ODxW0LAFSynTutFNriopWraWmhDmstIjFERZudDb8+P2+bX6iWDwCb09RkP6hN9AVQ0XbtCr++cydJIQAAm0ViiIp26FD4dcrjAACweUwlo6Lt32+bTq5etXUxO3daUjgyUuqWAQBQeUgMUfGGh62cwtqalclxrtQtAgCgMpEYoio4JzU0lLoVAABUtrzWGDrn/olz7qJz7rRz7r845zoK1TAAqFbEThRLcETp6dO25GZlpdQtQrnLd/PJ25Ke9t4/K+mypF/Jv0kAUPWIndh2S0vSt75lSeHNm9KFC9Kf/7kdCABkkldi6L3/mvd+LfbwXUn78m8SAFQ3YieK4dIlq+maaHVVOnu2NO1BZShkuZq/LemrBbwfANQCYie2xf374denp5lSRmYbbj5xzn1dUl/IU1/03v9R7DVflLQm6ctZ7vOWpLckaf/+/VtqLABUCmInSi3TcXbOcdQdMtswMfTefzrb8865n5f0OUmf8t77LPf5kqQvSdLx48czvg4AqgGxE6U2MCBdu5Z+vbfXSnsBYfLdlfxZSV+Q9NPe+4XCNAkAqhuxE8Vw5Ii0e3fytZYW6dlnS9MeVIZ8+wy/KalR0tvOqgq/673/xbxbBQDVjdiJbVdfL732mjQ1ZTuRd+2y0UIOAUA2eSWG3vsMJ9UC+VldtfpbTU2lbglQeMROFNPu3ekjh0AmrDJAWVletppbDx5I3kutrdJTT0k9PaVuGQAA1Y99SSgrH3xgJRaCpfhzc9L770uPH5e2XQAA1AISQ5SNhw/tJ1U0alX7AQDA9iIxRNlYyLI3M9tzAACgMEgMUTY6OjI/195evHYAAFCrSAxRNnbtsoKsqXbulIaGit8eAABqDbuSUVaef15qa5Nu3bKSNb29VqS1sbHULQMAoPqRGKKsOCcdOmQ/AACguEgMAQCoURMT0vy8zdR0dZW6NSgHJIYAANSY5WXp3Xel2dn4te5u6cUX7Sg91C42nwAAUGPOnk1OCiVpclK6dKk07UH5IDEEAKCGrK9L9+6FP3fnTnHbgvJDYggAQA3xPn7saKr19eK2BeWHxBAAgBpSX595o8mePcVtC8oPiSEAADXm6aelhobkazt3SseOlaY9KB/sPQIAoMa0t0uf/KQ0NhYvVzM4mJ4sovaQGAIAUIMaG6XDh0vdCpQbppIBAAAgicQQFWhtTZqZkZaWSt0SAACqC1PJqChXrkhXr1py6JzU1yd97GNU6gcAoBAYMUTFuHNHunjRkkLJ6nDduyedPl3adgEAUC1IDFExRkfDr9+9K62sFLUpAABUJRJDVIxMawq9JzEEAKAQSAxRMTJV6m9slJqbi9sWAACqEYkhKsbhw+HFV48dkyL8PxkAgLyxlxMVo6VFev1125U8M2PHNx04IPX2lrplAABUBxJDVJRdu6Tnnit1KwAAqE5MwAEAAEASiSEAAABiSAwBAAAgicQQAAAAMSSGAAAAkERiCAAAgBgSQwAAAEgiMQQAAEAMiSEAAAAkkRgCAAAghsQQAAAAkkgMAQAAEENiCAAAAEkkhgAAAIghMQQAAIAkEkMAAADEkBgCAABAEokhAAAAYkgMAQAAIKlAiaFz7vPOOe+c6y7E/QCgFhA7AZSbvBND59ygpB+TNJZ/cwCgNhA7AZSjQowY/jNJX5DkC3AvAKgVxE4AZSevxNA59zOS7njvPypQewCg6hE7AZSr+o1e4Jz7uqS+kKe+KOlXZVMhG3LOvSXpLUnav3//JpoIAJWH2AmgEjnvtzaL4Zx7RtL/J2khdmmfpLuSXvLe38/23uPHj/uTJ09u6XMB1Dbn3Cnv/fFSt2OriJ0ASiHX2LnhiGEm3vszknoTPnBU0nHv/eRW7wkA1Y7YCaCcUccQAAAAkvIYMUzlvR8u1L0AoFYQOwGUE0YMAQAAIInEEAAAADEkhgAAAJBEYggAAIAYEkMAAABIIjEEAABADIkhAAAAJJEYAgAAIIbEEAAAAJJIDAEAABBDYggAAABJJIYAAACIITEEAACAJBJDAAAAxJAYAgAAQBKJIQAAAGJIDAEAACCJxBAAAAAxJIYAAACQRGIIAACAGBJDAAAASCIxBAAAQIzz3hf/Q52bkHSzSB/XLWmySJ9VTNX6vaTq/W7V+r2k4n63Ie99T5E+q6wQOwuiWr+XVL3fje9VGDnFzpIkhsXknDvpvT9e6nYUWrV+L6l6v1u1fi+pur9brarW/6bV+r2k6v1ufK/iYioZAAAAkkgMAQAAEFMLieGXSt2AbVKt30uq3u9Wrd9Lqu7vVquq9b9ptX4vqXq/G9+riKp+jSEAAAByUwsjhgAAAMhBTSWGzrnPO+e8c6671G0pBOfcP3HOXXTOnXbO/RfnXEep25QP59xnnXOXnHNXnXP/S6nbUyjOuUHn3Decc+edc+ecc79U6jYVknOuzjn3fefcn5a6LdgexM7yRuysTOUaO2smMXTODUr6MUljpW5LAb0t6Wnv/bOSLkv6lRK3Z8ucc3WSfkvST0h6UtJfc849WdpWFcyapM9775+U9Iqk/7GKvpsk/ZKkC6VuBLYHsbO8ETsrWlnGzppJDCX9M0lfkFQ1iyq991/z3q/FHr4raV8p25OnlyRd9d5f996vSPo9ST9T4jYVhPf+nvf+w9if52SBYKC0rSoM59w+ST8l6XdK3RZsG2JneSN2VqByjp01kRg6535G0h3v/Uelbss2+tuSvlrqRuRhQNKthMe3VSUBIJFzbljS85LeK21LCuY3ZElDtNQNQeEROysCsbMylW3srC91AwrFOfd1SX0hT31R0q/KpkIqTrbv5b3/o9hrvigbcv9yMduGzXHOtUj6z5J+2Xs/W+r25Ms59zlJ4977U865N0rdHmwNsZPYWe6IncVVNYmh9/7TYdedc89IOiDpI+ecZFMGHzrnXvLe3y9iE7ck0/cKOOd+XtLnJH3KV3btoTuSBhMe74tdqwrOuQZZYPuy9/4PS92eAjkh6aedcz8paaekNufcv/fe//UStwubQOwkdpYzYmfx1VwdQ+fcqKTj3vuKP5DbOfdZSf9U0o967ydK3Z58OOfqZYvAPyULah9I+jnv/bmSNqwAnP1W/V1J0977Xy51e7ZDrNf7D733nyt1W7A9iJ3lidhZ2coxdtbEGsMq9puSWiW97Zz7gXPut0vdoK2KLQT/+5L+TLbA+D9WQ2CLOSHpb0h6M/bf6QexniKA0iB2VgZiZwnU3IghAAAAwjFiCAAAAEkkhgAAAIghMQQAAIAkEkMAAADEkBgCAABAEokhAAAAYkgMAQAAIInEEAAAADH/P499us20qGEfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 792x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# cf. original normal randoms in cell [2]\n",
    "plt.figure(figsize=(11,5)); \n",
    "plt.subplot(1,2,1).set_aspect('equal'); plt.xlim(-5,5); plt.ylim(-5,5); plt.title('whitened')\n",
    "plt.scatter(B[0,:],B[1,:], marker='o',color='b', s=50, alpha=0.3, edgecolor='none');\n",
    "plt.subplot(1,2,2).set_aspect('equal'); plt.xlim(-5,5); plt.ylim(-5,5); plt.title('original')\n",
    "plt.scatter(-N[0,:],N[1,:], marker='o',color='b', s=50, alpha=0.3, edgecolor='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Whitening using sklearn package (same result)\n",
    "- Set \"whiten\" to \"True\" when initializing your PCA instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'PC2')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAAFACAYAAAAruW7uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEdpJREFUeJzt3f2PXOV5xvHr8q6NN2BitzZYsHZNA1VKgRJl40ZFVRJMUkhIEFIV0TaRkqi1UjWSkahQgD+gPyAlREqkyApVowbJaUSiRJQ2MRAq+gPBa2JwDYQAIXh5CSbY2IJV7LXv/rCz9GE8u7PzcuY5L9+PZGl3ZnbnPmvttfdznuc8xxEhAMC8FbkLAIAyIRQBIEEoAkCCUASABKEIAAlCEQAShCIAJAhFAEgQigCQGM9dQC/Wr18fW7ZsyV0GgArau3fvaxGxodvrKhWKW7Zs0fT0dO4yAFSQ7V8v53UMnwEgQSgCQIJQBIAEoQgACUIRABKEIgAkCEUASBCKAJAgFAEgQSgCQIJQBIAEoQgACUIRABKEIgAkCEUASBCKAJAgFAEgQSgCQIJQBIAEoQgACUIRABKEIgAksoei7THbP7d9T+5aACB7KEraIenJ3EUAgJQ5FG1PSvqEpG/lrAMAFuTuFO+QdLOkU4u9wPZ229O2pw8dOjS6ygA0UrZQtH2tpFcjYu9Sr4uInRExFRFTGzZsGFF1AJoqZ6d4haRP2X5e0i5JV9r+TsZ6ACBfKEbELRExGRFbJN0g6YGI+EyuegBAyn9OEQBKZTx3AZIUEQ9KejBzGQBApwgAKUIRABKEIgAkCEUASBCKAJAgFAEgQSgCQIJQBIAEoQgACUIRABKEIgAkCEUASBCKAJAgFAEgQSgCQIJQBIAEoQgACUIRABKEIgAkCEWU2ktHZvXIr17XS0dmc5eChijFjauATl46Mquv3f9LzZ08pfGxFdqx7SKdt3Yid1moOTpFlNbM4VnNnTylyXXv0tzJU5o5TLeI4hGKKK3JdRMaH1uhmcNvaXxshSbX0SWieAyfUVrnrZ3Qjm0XaebwrCbXTTB0xkgQiii189YShhgths8AkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhGVxu0KMGxsHYbK4nYFKAKdIiqL2xWgCIQiKovbFaAIDJ9RWdyuAEUgFFFp3K4Aw8bwGQAS2ULR9ibbP7X9hO0DtnfkqgUAFuQcPs9JuikiHrW9RtJe27sj4omMNQFouGydYkS8HBGPtj4+JulJSefnqgcApJKcU7S9RdL7JP0sbyUAmi57KNo+S9Ldkm6MiKMdnt9ue9r29KFDh0ZfIIBGyRqKtldqPhDviojvd3pNROyMiKmImNqwYcNoCwTQODlnny3pTklPRsRXctUBAKmcneIVkj4r6Urb+1r/Pp6xHgDItyQnIv5HknO9PwB0kn2iBQDKhFAEgAShiNJgF22UAbvkoBTYRRtlQaeIUmAXbZQFoYhSYBdtlAXDZ5QCu2ijLAhFlAa7aKMMGD4DQIJQBIAEoYjGYl0kOuGcIhqJdZFYDJ0iGol1kVgMoYhGYl0kFsPwGY1U9nWRLx2ZLW1tdUcoorG/gGVdF8n5zrwIxYbjF7B80vOdM4ff0szhWf5PRohzig3HhMNgiljWw/nOvOgUG45fwP4V1WWX/Xxn3RGKDccvYP+KHOaW9XxnExCK4BewT3TZ9UQoAn2iy64nQhFYpk5Ll+iy64dQBJaBpUvNwZIc1E4Ry2RYutQcdIqolaI6OiZVmoNQRK0UtUyGSZXmIBRRK0V2dEyqNAOhiFqho8OgCEXUDh0dBsHsM/rGPU5QR3SK6Avr9lBXdIroS451e3SmGAU6RfSl11neQXf3LmNn2tQdy+uOUERf2md5JemRX73+9sdpWAwj0Ja7/nBUQVXGkMZwEIro28IsbxoQv5s7JUtaNb7i7bAYxoLq5XSmgwZVL4HKLQPqi1DEwNKA2HfwiKTQ5ZvWvR0Ww1hQvdj6wzTIBgmqXgOVy/7qi1DEwNKAWLN6XJbeERbDWlDdvv6wPchu+MCmvoOq10BlkXh9EYroWfsws9P5xVHsO9geZCdORt9B1U/nxyLxeiIU0ZPFhpntATGKsOgUZP0GFZ0fFhCK6EmZJhiGHWR0fpAIRfRolBMM7cN0bgeAUegairbPlrQhIp5te/yyiHi8sMpQSqMaZnaaRNm15yDrAlG4JS/zs/1pSU9Jutv2AdsfSJ7+1yILQ3mdt3ZCWy/4vUJDqf0ywsdn3uB2ABiJbtc+3yrp/RFxuaTPS/o329e3nnOhlSGbMlxj3D5Mv2zy3QMN28twTKiGbsPnsYh4WZIi4hHbH5F0j+1NkqLw6jByZbl8rdMw/dyzV/c1bF/OMXEdMxZ06xSP2X7PwietgPywpOsk/cmgb277atu/sP2M7S8P+v0wuDLdta59mN7vsL3bMS2E5nf3vKCv3f/LZXWTdJ711a1T/Ae1DZMj4pjtqyV9epA3tj0m6RuSPippRtIe2z+KiCcG+b4YTJGzy7m6sW7H1Osyo7J00yhGt1B8U9K5kp5pe3yrpIcHfO+tkp6JiOckyfYuzXeghGJGRc0u5wySbsfU6x+CMq3VxPB1C8U7JN3S4fGjrec+OcB7ny/pYPL5jKQ/a3+R7e2StkvS5s2bB3g7LNcoLskbdZAsdUy9/iFgM4h66xaK50bE/vYHI2K/7S2FVHT6e+2UtFOSpqammNypqH6CZJTD7V7+ECw3RIusn4mh4nQLxbVLPDfo/8SLkjYln0+2HkMN9dqNlf28XbcQLbL+sv9sqq7b7PO07b9vf9D230naO+B775F0ke0LbK+SdIOkHw34PVEi7TO0vcwel2kWvB/91L/cGe2q/2zKrluneKOkH9j+W/1/CE5JWiXp+kW/ahkiYs72lyT9WNKYpH+JiAODfE+Ux6DdTNXP2/VzD5vl/ryq/rMpuyVDMSJ+I+nPW4u2L2k9/B8R8cAw3jwi7pV07zC+F8pl0ImVMmzlNch5u17r7+XnVYafTZ0tGYq2V0v6oqQLJe2XdGdEzI2iMFTbsG5BkOsXfhjn7Xqpv9efF7sDFafb8Pnbkk5IekjSNZL+WPNDamBJRXQzo5xxHfUSIrq/8ugWihdHxKWSZPtOSY8UXxLqYpjdzKhnXHOct6P7K4duoXhi4YPWxEjB5QCdzRye1dHZ4zrzjJU6Onuczg2F6RaKf2r7aOtjS5pofW5JERFnF1od0LJyzHrqlWM6eSo0tsJaOVb8H2g6t2bqNvs8NqpC0Cy9nh88cTL03o1rdOaqlXrz+AmdOFnfi5u4WiUv7tGCkevn/ODkugmdPbFKcydP6eyJVQOf4ytr8HC1Sn6EIkaun5ndYZ7jK3Pw5N44A4Rio5V1f8PFDOscX5mDh6tV8iMUG6rM+xsWrczBk/tnA0KxsXJ3SzlndssePMx650UoNlSZu6VRIHiwGEKxocreLQG5EIoNRrcEnK7bJrMA0CiEIgAkCEU0CjexRzecU0RjlPlKFpQHnSIao/2GT/sOHqFrxGnoFNEY6drM43OndO/+l3XG+Iolu8aybhyB4hCKaIx0bebTvzmm3U+8onM2nKU3Zk90vKKH4XYzMXxGLSx3AuW8tROaXDehh5/7rZ7/7Vu6/6lXdXzulFaO+bSv5/7KzUSniI6qNGzstaObOTyrM8ZXaNt7z9Gzh97UB//w97Vrz8HTvr7pl0I2FaGI01Rt2Njr5hYLYffG7AltfPdqrTtzVcev51LIZiIUcZrcO+j0anLdhH43Nz+bvGb1+LLumZyGnST999OHOnaEXArZPIQiTlPFYeP8baxCy72dVXvY0RFiAaGI01Rt2DhzeFarxlfo8k3r+u5s6QixgFBER1UKiSp2tigvQhGVV7XOFuVGKKIWqtTZotxYvA0ACUIRABKEIgAkCEUASBCKAJAgFFFq3D4Ao8aSHJRW1TamQD3QKaK0hrGfIZ0mekWniNIa9PI9Ok30g1BEafVy+V6nTXGrtgUayoFQRKkt5/K9xTpCNopAPwhFVN5iHSEbRaAfhCIqb6mOkI0i0CtCEZVHR4hhIhTRtzLd8Y+OEMOSJRRt3y7pk5KOS3pW0ucj4kiOWtAflrugrnIt3t4t6ZKIuEzS05JuyVQH+sSN4lFXWUIxIn4SEXOtTx+WNJmjDvSP5S6oqzKcU/yCpO8u9qTt7ZK2S9LmzZtHVRO6YHIDdVVYKNq+T9LGDk/dFhE/bL3mNklzku5a7PtExE5JOyVpamoqCigVfWJyA3VUWChGxFVLPW/7c5KulbQtIgg7AKWQa/b5akk3S/pQRLyVowYA6CTX7PPXJa2RtNv2PtvfzFQHALxDlk4xIi7M8b4A0A2bzAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACQIRQBIEIoAkCAUASBBKAJAglAEgAShCAAJQhEAEoQiACSyhqLtm2yH7fU56wCABdlC0fYmSR+T9EKuGgCgXc5O8auSbpYUGWsAgHfIEoq2r5P0YkQ8luP9AWAx40V9Y9v3SdrY4anbJN2q+aHzcr7PdknbJWnz5s1Dqw8AOnHEaEevti+VdL+kt1oPTUp6SdLWiHhlqa+dmpqK6enpgisEUEe290bEVLfXFdYpLiYi9ks6Z+Fz289LmoqI10ZdCwC0Y50iACRG3im2i4gtuWsAgAV0igCQIBQBIEEoAkCCUASABKEIAAlCEQAShCIAJAhFAEgQigCQIBQBIEEoAkCCUASABKEIAAlCEQAShCIAJAhFAEgQigCQIBQBIEEoAkCCUASABKEIAAlCEQASjojcNSyb7UOSfj2it1sv6bURvdco1fW4pPoeW12PSxrtsf1BRGzo9qJKheIo2Z6OiKncdQxbXY9Lqu+x1fW4pHIeG8NnAEgQigCQIBQXtzN3AQWp63FJ9T22uh6XVMJj45wiACToFAEgQSgCQIJQXAbbN9kO2+tz1zIMtm+3/ZTtx23/wPba3DUNwvbVtn9h+xnbX85dz7DY3mT7p7afsH3A9o7cNQ2T7THbP7d9T+5aUoRiF7Y3SfqYpBdy1zJEuyVdEhGXSXpa0i2Z6+mb7TFJ35B0jaSLJf217YvzVjU0c5JuioiLJX1Q0j/W6NgkaYekJ3MX0Y5Q7O6rkm6WVJsZqYj4SUTMtT59WNJkznoGtFXSMxHxXEQcl7RL0nWZaxqKiHg5Ih5tfXxM8wFyft6qhsP2pKRPSPpW7lraEYpLsH2dpBcj4rHctRToC5L+M3cRAzhf0sHk8xnVJDhStrdIep+kn+WtZGju0HyzcSp3Ie3GcxeQm+37JG3s8NRtkm7V/NC5cpY6roj4Yes1t2l+iHbXKGtDb2yfJeluSTdGxNHc9QzK9rWSXo2IvbY/nLuedo0PxYi4qtPjti+VdIGkx2xL80PMR21vjYhXRlhiXxY7rgW2PyfpWknbotqLVV+UtCn5fLL1WC3YXqn5QLwrIr6fu54huULSp2x/XNJqSWfb/k5EfCZzXZJYvL1stp+XNBURld+txPbVkr4i6UMRcSh3PYOwPa75yaJtmg/DPZL+JiIOZC1sCDz/1/jbkl6PiBtz11OEVqf4TxFxbe5aFnBOsZm+LmmNpN2299n+Zu6C+tWaMPqSpB9rfiLi3+sQiC1XSPqspCtb/0/7Wt0VCkSnCAAJOkUASBCKAJAgFAEgQSgCQIJQBIAEoYhKsH2ytSTlf21/z/a7Wo9vtL3L9rO299q+1/YftZ77L9tHyrYLC8qNUERVzEbE5RFxiaTjkr7YWtz8A0kPRsR7IuL9mt/x59zW19yu+XV+wLIRiqiihyRdKOkjkk5ExNuLzyPisYh4qPXx/ZKO5SkRVUUoolJal/VdI2m/pEsk7c1bEeqGUERVTNjeJ2la8xv+3pm5HtRU43fJQWXMRsTl6QO2D0j6q0z1oKboFFFlD0g6w/b2hQdsX2b7LzLWhIojFFFZrX0gr5d0VWtJzgFJ/yzpFUmy/ZCk70naZnvG9l/mqxZVwS45AJCgUwSABKEIAAlCEQAShCIAJAhFAEgQigCQIBQBIPF/32gXpoiAkkUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sklearn\n",
    "# Whiten the data and choose first 4 components\n",
    "from sklearn.decomposition import PCA\n",
    "pca = PCA(n_components=X[:,0].size, whiten = True).fit_transform(X.T) # The only difference is Whiten = True\n",
    "\n",
    "#plot the first two principal components\n",
    "plt.figure(figsize=(5,5)); plt.subplot(111,aspect='equal');plt.xlim(-5,5); plt.ylim(-5,5);\n",
    "plt.plot(pca[:,0],pca[:,1], '.', alpha=0.5); plt.xlabel('PC1'); plt.ylabel('PC2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scree Plot\n",
    "\n",
    "- The eigenvalue spectrum\n",
    "\n",
    ">$ \\big\\{ \\lambda_1, \\lambda_2, \\dots, \\lambda_N \\big\\}$\n",
    "\n",
    "- How many important directions?\n",
    "\n",
    "> Keep $K =\\,?$ principal components\n",
    "\n",
    "- Explained variance \n",
    "\n",
    "> Cf. $\\mathbb{Var}[X\\pm{}Y] = \\mathbb{Var}[X]+\\mathbb{Var}[Y]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate 10-D vectors: scale, rotate\n",
    "np.random.seed(1)\n",
    "Z = gaussian.rvs(0,1,(10,1000))\n",
    "if True: # scale them here\n",
    "    for i in range(Z[:,0].size): \n",
    "        Z[i,:] *= np.sqrt(i)\n",
    "    Z[:4,:] *= 1e-7\n",
    "# quick-n-dirty random rotation\n",
    "M = np.random.randn(Z[:,0].size,Z[:,0].size)\n",
    "Q,_ = np.linalg.qr(M) # QR decomposition\n",
    "Y = Q.dot(Z) # random rotation\n",
    "print (Y.shape)\n",
    "np.savetxt(\"temp.csv\", Y.T, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove all previous variables from memory\n",
    "del Y, M, Q, Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pandas dataframe - table data structure\n",
    "df = pd.read_csv('temp.csv',header=None)\n",
    "df[:3] # slice like arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# re-load the Y matrix and convert dataframe to matrix\n",
    "Y = pd.read_csv('temp.csv',header=None).values.T\n",
    "Y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(111,aspect='equal')\n",
    "plt.plot(Y[0,:],Y[1,:], '.', alpha=0.2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import decomposition\n",
    "pca = decomposition.PCA(n_components=Y.shape[0])\n",
    "B = pca.fit_transform(Y).T\n",
    "E, L = pca.components_.T, pca.explained_variance_\n",
    "\n",
    "plt.subplot(111,aspect='equal')\n",
    "plt.plot(B[0,:],B[5,:], '.', alpha=0.2);\n",
    "# something is wrong - fix the bug above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(211); plt.plot(L,'o-'); plt.ylabel('Eigenvalues');\n",
    "plt.subplot(212); cl = np.cumsum(L); plt.ylabel('Total Variance');\n",
    "plt.plot(cl,'o-r'); plt.ylim(0,None);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot orig and new covariance matrices (estimate w/o norm)\n",
    "plt.figure(figsize=(5,2.5))\n",
    "plt.subplot(121); plt.imshow(Y.dot(Y.T),interpolation='none');\n",
    "plt.subplot(122); plt.imshow(B.dot(B.T),interpolation='none');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = B[:4,:]\n",
    "# plot orig and new covariance matrices (estimate w/o norm)\n",
    "plt.figure(figsize=(5,2.5))\n",
    "plt.subplot(121); plt.imshow(B.dot(B.T),interpolation='none');\n",
    "plt.subplot(122); plt.imshow(A.dot(A.T),interpolation='none');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inverse of the Covariance Matrix\n",
    "\n",
    "- Appears in the multivariate normal distribution!\n",
    "\n",
    ">$\\displaystyle{\\cal{}N}(x;\\mu,C) \\propto \\exp\\left[-\\frac{1}{2}(x\\!-\\!\\mu)^T\\,C^{-1} (x\\!-\\!\\mu)\\right]$\n",
    "\n",
    "- Inverse of the diagonal eigenvalue matrix\n",
    "\n",
    ">$\\displaystyle \\Lambda^{-1} =  \\left( \\begin{array}{ccc}\n",
    "\\frac{1}{\\lambda_1} &  & \\cdots & 0\\\\\n",
    " & \\frac{1}{\\lambda_2} &   & \\vdots\\\\\n",
    "\\vdots &  & \\ddots &  \\\\\n",
    "0 & \\cdots &  & \\frac{1}{\\lambda_N} \\\\\n",
    "\\end{array} \\right)$\n",
    "\n",
    "- Inverse of the covariance matrix\n",
    "\n",
    ">$\\displaystyle C^{-1} = E\\ \\Lambda^{-1} E^T$\n",
    "\n",
    "- Also see pseudoinverse with small eigenvalues "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting Lines\n",
    "\n",
    "- What if $x$ and $y$ are both noisy? \n",
    "\n",
    "> For example, $\\big\\{(x_i,y_i)\\big\\}$ measurements have the same uncertainties. \n",
    "> The relevant residuals are perpendicular to the line.\n",
    "> Minimizing RMS of residuals is related to maximizing the sample variance along line!\n",
    "\n",
    "- Sounds like the PCA problem?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting Planes\n",
    "\n",
    "- Similarly, fitting a $K$-dimensional hyperplane in $N$ dimensions\n",
    "\n",
    "> Minimizing sum of square lengths of the residual vectors\n",
    "\n",
    ">$\\displaystyle \\ \\ \\ \\ \\ \\min \\sum_i r_i^2 \\ \\ \\ \\ \\ $  where $\\ \\ \\ r_i = x_i - (a\\,a^T)x_i$, \n",
    "\n",
    "> yields \n",
    "\n",
    ">$\\displaystyle \\ \\ \\ \\ \\ \\max \\sum_i \\left(a^Tx_i\\right)^2 \\ \\ \\ \\ \\ $ \n",
    "\n",
    "> cf. sample variance along $a$, if data already center\n",
    "\n",
    "- Essentially same as the PCA problem!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate 2D (column) vectors\n",
    "np.random.seed(seed=42);\n",
    "N = gaussian.rvs(0,1,(2,50)); \n",
    "N[0,:] *= 2 \n",
    "f = +pi/6   # rotate by 30 deg\n",
    "R = np.array([[cos(f), -sin(f)],\n",
    "           [sin(f),  cos(f)]]) \n",
    "X = R.dot(N)\n",
    "\n",
    "plt.figure(figsize=(5,5)); plt.subplot(111,aspect='equal');\n",
    "plt.plot(X[0,:],X[1,:],'x',alpha=0.6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# project on 1st pricipal component\n",
    "X -= np.mean(X, axis=1).reshape(X[:,1].size,1)\n",
    "E,_,_ = np.linalg.svd(X) # only eigenvectors\n",
    "F = E[:,:1]           # truncated basis: only PC1\n",
    "P = F.dot(F.T).dot(X) # projection\n",
    "R = X - P;            # residuals\n",
    "\n",
    "plt.figure(figsize=(5,5)); plt.subplot(111,aspect='equal');\n",
    "plt.plot(X[0,:],X[1,:],'xb',alpha=0.6);\n",
    "plt.plot(P[0,:],P[1,:],'-or',alpha=0.4);\n",
    "plt.quiver(P[0,:],P[1,:],R[0,:],R[1,:], alpha=0.2,\n",
    "    angles='xy',scale_units='xy',scale=1);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
